ReadNCDF.cpp

Go to the documentation of this file.

/**
 * MOAB, a Mesh-Oriented datABase, is a software component for creating,
 * storing and accessing finite element mesh data.
 *
 * Copyright 2004 Sandia Corporation.  Under the terms of Contract
 * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
 * retains certain rights in this software.
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 *
 */

#ifdef WIN32
#pragma warning( disable : 4786 )
#endif

#include "ReadNCDF.hpp"
#include "netcdf.h"

#include <algorithm>
#include <string>
#include <cassert>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <sstream>
#include <iostream>
#include <map>

#include "moab/CN.hpp"
#include "moab/Range.hpp"
#include "moab/Interface.hpp"
#include "ExoIIUtil.hpp"
#include "MBTagConventions.hpp"
#include "Internals.hpp"
#include "moab/ReadUtilIface.hpp"
#include "exodus_order.h"
#include "moab/FileOptions.hpp"
#include "moab/AdaptiveKDTree.hpp"
#include "moab/CartVect.hpp"

namespace moab
{

#define INS_ID( stringvar, prefix, id ) sprintf( stringvar, prefix, id )

#define GET_DIM( ncdim, name, val )                                                   \
    {                                                                                 \
        int gdfail = nc_inq_dimid( ncFile, name, &( ncdim ) );                        \
        if( NC_NOERR == gdfail )                                                      \
        {                                                                             \
            size_t tmp_val;                                                           \
            gdfail = nc_inq_dimlen( ncFile, ncdim, &tmp_val );                        \
            if( NC_NOERR != gdfail )                                                  \
            {                                                                         \
                MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Couldn't get dimension length" ); \
            }                                                                         \
            else                                                                      \
                ( val ) = tmp_val;                                                    \
        }                                                                             \
        else                                                                          \
            ( val ) = 0;                                                              \
    }

#define GET_DIMB( ncdim, name, varname, id, val ) \
    INS_ID( name, varname, id );                  \
    GET_DIM( ncdim, name, val );

#define GET_VAR( name, id, dims )                                                               \
    {                                                                                           \
        ( id )     = -1;                                                                        \
        int gvfail = nc_inq_varid( ncFile, name, &( id ) );                                     \
        if( NC_NOERR == gvfail )                                                                \
        {                                                                                       \
            int ndims;                                                                          \
            gvfail = nc_inq_varndims( ncFile, id, &ndims );                                     \
            if( NC_NOERR == gvfail )                                                            \
            {                                                                                   \
                ( dims ).resize( ndims );                                                       \
                gvfail = nc_inq_vardimid( ncFile, id, &( dims )[0] );                           \
                if( NC_NOERR != gvfail )                                                        \
                {                                                                               \
                    MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Couldn't get variable dimension IDs" ); \
                }                                                                               \
            }                                                                                   \
        }                                                                                       \
    }

#define GET_1D_INT_VAR( name, id, vals )                                                      \
    {                                                                                         \
        GET_VAR( name, id, vals );                                                            \
        if( -1 != ( id ) )                                                                    \
        {                                                                                     \
            size_t ntmp;                                                                      \
            int ivfail = nc_inq_dimlen( ncFile, ( vals )[0], &ntmp );                         \
            if( NC_NOERR != ivfail )                                                          \
            {                                                                                 \
                MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Couldn't get dimension length" );         \
            }                                                                                 \
            ( vals ).resize( ntmp );                                                          \
            size_t ntmp1 = 0;                                                                 \
            ivfail       = nc_get_vara_int( ncFile, id, &ntmp1, &ntmp, &( vals )[0] );        \
            if( NC_NOERR != ivfail )                                                          \
            {                                                                                 \
                MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting variable " << ( name ) ); \
            }                                                                                 \
        }                                                                                     \
    }

#define GET_1D_DBL_VAR( name, id, vals )                                                      \
    {                                                                                         \
        std::vector< int > dum_dims;                                                          \
        GET_VAR( name, id, dum_dims );                                                        \
        if( -1 != ( id ) )                                                                    \
        {                                                                                     \
            size_t ntmp;                                                                      \
            int dvfail = nc_inq_dimlen( ncFile, dum_dims[0], &ntmp );                         \
            if( NC_NOERR != dvfail )                                                          \
            {                                                                                 \
                MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Couldn't get dimension length" );         \
            }                                                                                 \
            ( vals ).resize( ntmp );                                                          \
            size_t ntmp1 = 0;                                                                 \
            dvfail       = nc_get_vara_double( ncFile, id, &ntmp1, &ntmp, &( vals )[0] );     \
            if( NC_NOERR != dvfail )                                                          \
            {                                                                                 \
                MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting variable " << ( name ) ); \
            }                                                                                 \
        }                                                                                     \
    }

ReaderIface* ReadNCDF::factory( Interface* iface )
{
    return new ReadNCDF( iface );
}

ReadNCDF::ReadNCDF( Interface* impl ) : mdbImpl( impl ), max_line_length( -1 ), max_str_length( -1 )
{
    assert( impl != NULL );
    reset();

    impl->query_interface( readMeshIface );

    // Initialize in case tag_get_handle fails below
    mMaterialSetTag  = 0;
    mDirichletSetTag = 0;
    mNeumannSetTag   = 0;
    mHasMidNodesTag  = 0;
    mDistFactorTag   = 0;
    mQaRecordTag     = 0;
    mGlobalIdTag     = 0;

    //! Get and cache predefined tag handles
    ErrorCode result;
    const int negone = -1;
    result           = impl->tag_get_handle( MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mMaterialSetTag,
                                             MB_TAG_SPARSE | MB_TAG_CREAT, &negone );
    assert( MB_SUCCESS == result );
    result = impl->tag_get_handle( DIRICHLET_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mDirichletSetTag,
                                   MB_TAG_SPARSE | MB_TAG_CREAT, &negone );
    assert( MB_SUCCESS == result );
    result = impl->tag_get_handle( NEUMANN_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mNeumannSetTag,
                                   MB_TAG_SPARSE | MB_TAG_CREAT, &negone );
    assert( MB_SUCCESS == result );
    const int mids[] = { -1, -1, -1, -1 };
    result           = impl->tag_get_handle( HAS_MID_NODES_TAG_NAME, 4, MB_TYPE_INTEGER, mHasMidNodesTag,
                                             MB_TAG_SPARSE | MB_TAG_CREAT, mids );
    assert( MB_SUCCESS == result );
    result = impl->tag_get_handle( "distFactor", 0, MB_TYPE_DOUBLE, mDistFactorTag,
                                   MB_TAG_SPARSE | MB_TAG_VARLEN | MB_TAG_CREAT );
    assert( MB_SUCCESS == result );
    result = impl->tag_get_handle( "qaRecord", 0, MB_TYPE_OPAQUE, mQaRecordTag,
                                   MB_TAG_SPARSE | MB_TAG_VARLEN | MB_TAG_CREAT );
    assert( MB_SUCCESS == result );
    mGlobalIdTag = impl->globalId_tag();

    assert( MB_SUCCESS == result );
#ifdef NDEBUG
    if( MB_SUCCESS == result )
    {
    };  // Line to avoid compiler warning about unused variable
#endif
    ncFile = 0;
}

void ReadNCDF::reset()
{
    mCurrentMeshHandle = 0;
    vertexOffset       = 0;

    numberNodes_loading         = 0;
    numberElements_loading      = 0;
    numberElementBlocks_loading = 0;
    numberFaceBlocks_loading    = 0;
    numberNodeSets_loading      = 0;
    numberSideSets_loading      = 0;
    numberDimensions_loading    = 0;

    if( !blocksLoading.empty() ) blocksLoading.clear();

    if( !nodesInLoadedBlocks.empty() ) nodesInLoadedBlocks.clear();

    if( !faceBlocksLoading.empty() ) faceBlocksLoading.clear();
}

ReadNCDF::~ReadNCDF()
{
    mdbImpl->release_interface( readMeshIface );
}

ErrorCode ReadNCDF::read_tag_values( const char* file_name,
                                     const char* tag_name,
                                     const FileOptions& /* opts */,
                                     std::vector< int >& id_array,
                                     const SubsetList* subset_list )
{
    if( subset_list )
    {
        MB_SET_ERR( MB_UNSUPPORTED_OPERATION, "ExodusII reader supports subset read only by material ID" );
    }

    // Open netcdf/exodus file
    int fail = nc_open( file_name, 0, &ncFile );
    if( NC_NOWRITE != fail )
    {
        MB_SET_ERR( MB_FILE_DOES_NOT_EXIST, "ReadNCDF:: problem opening Netcdf/Exodus II file " << file_name );
    }

    // 1. Read the header
    ErrorCode rval = read_exodus_header();
    if( MB_FAILURE == rval ) return rval;

    int count = 0;
    const char* prop;
    const char* blocks   = "eb_prop1";
    const char* nodesets = "ns_prop1";
    const char* sidesets = "ss_prop1";

    if( !strcmp( tag_name, MATERIAL_SET_TAG_NAME ) )
    {
        count = numberElementBlocks_loading;
        prop  = blocks;
    }
    else if( !strcmp( tag_name, DIRICHLET_SET_TAG_NAME ) )
    {
        count = numberNodeSets_loading;
        prop  = nodesets;
    }
    else if( !strcmp( tag_name, NEUMANN_SET_TAG_NAME ) )
    {
        count = numberSideSets_loading;
        prop  = sidesets;
    }
    else
    {
        ncFile = 0;
        return MB_TAG_NOT_FOUND;
    }

    if( count )
    {
        int nc_var = -1;
        GET_1D_INT_VAR( prop, nc_var, id_array );
        if( !nc_var )
        {
            MB_SET_ERR( MB_FAILURE, "Problem getting prop variable" );
        }
    }

    // Close the file
    fail = nc_close( ncFile );
    if( NC_NOERR != fail )
    {
        MB_SET_ERR( MB_FAILURE, "Trouble closing file" );
    }

    ncFile = 0;
    return MB_SUCCESS;
}

ErrorCode ReadNCDF::load_file( const char* exodus_file_name,
                               const EntityHandle* file_set,
                               const FileOptions& opts,
                               const ReaderIface::SubsetList* subset_list,
                               const Tag* file_id_tag )
{
    ErrorCode status;
    int fail;

    int num_blocks            = 0;
    const int* blocks_to_load = 0;
    if( subset_list )
    {
        if( subset_list->tag_list_length > 1 || !strcmp( subset_list->tag_list[0].tag_name, MATERIAL_SET_TAG_NAME ) )
        {
            MB_SET_ERR( MB_UNSUPPORTED_OPERATION, "ExodusII reader supports subset read only by material ID" );
        }
        if( subset_list->num_parts )
        {
            MB_SET_ERR( MB_UNSUPPORTED_OPERATION, "ExodusII reader does not support mesh partitioning" );
        }
        blocks_to_load = subset_list->tag_list[0].tag_values;
        num_blocks     = subset_list->tag_list[0].num_tag_values;
    }

    // This function directs the reading of an exoii file, but doesn't do any of
    // the actual work

    // See if opts has tdata.
    ErrorCode rval;
    std::string s;
    rval = opts.get_str_option( "tdata", s );
    if( MB_SUCCESS == rval && !s.empty() )
        return update( exodus_file_name, opts, num_blocks, blocks_to_load, *file_set );

    reset();

    // 0. Open the file.

    // open netcdf/exodus file
    fail = nc_open( exodus_file_name, 0, &ncFile );
    if( NC_NOERR != fail )
    {
        MB_SET_ERR( MB_FILE_DOES_NOT_EXIST, "ReadNCDF:: problem opening Netcdf/Exodus II file " << exodus_file_name );
    }

    // 1. Read the header
    status = read_exodus_header();
    if( MB_FAILURE == status ) return status;

    status = mdbImpl->get_entities_by_handle( 0, initRange );
    if( MB_FAILURE == status ) return status;

    // 2. Read the nodes unless they've already been read before
    status = read_nodes( file_id_tag );
    if( MB_FAILURE == status ) return status;

    // 3.
    // extra for polyhedra blocks
    if( numberFaceBlocks_loading > 0 )
    {
        status = read_face_blocks_headers();
        if( MB_FAILURE == status ) return status;
    }

    status = read_block_headers( blocks_to_load, num_blocks );
    if( MB_FAILURE == status ) return status;

    // 4. Read elements (might not read them, depending on active blocks)
    if( numberFaceBlocks_loading > 0 )
    {
        status = read_polyhedra_faces();
        if( MB_FAILURE == status ) return status;
    }
    status = read_elements( file_id_tag );
    if( MB_FAILURE == status ) return status;

    // 5. Read global ids
    status = read_global_ids();
    if( MB_FAILURE == status ) return status;

    // 6. Read nodesets
    status = read_nodesets();
    if( MB_FAILURE == status ) return status;

    // 7. Read sidesets
    status = read_sidesets();
    if( MB_FAILURE == status ) return status;

    // 8. Read qa records
    if( file_set )
    {
        status = read_qa_records( *file_set );
        if( MB_FAILURE == status ) return status;
    }

    // What about properties???

    // Close the file
    fail = nc_close( ncFile );
    if( NC_NOERR != fail )
    {
        MB_SET_ERR( MB_FAILURE, "Trouble closing file" );
    }

    ncFile = 0;
    return MB_SUCCESS;
}

ErrorCode ReadNCDF::read_exodus_header()
{
    CPU_WORD_SIZE = sizeof( double );  // With ExodusII version 2, all floats
    IO_WORD_SIZE  = sizeof( double );  // should be changed to doubles

    // NetCDF doesn't check its own limits on file read, so check
    // them here so it doesn't corrupt memory any more than absolutely
    // necessary.
    int ndims;
    int fail = nc_inq_ndims( ncFile, &ndims );
    if( NC_NOERR != fail || ndims > NC_MAX_DIMS )
    {
        MB_SET_ERR( MB_FAILURE, "ReadNCDF: File contains " << ndims << " dims but NetCDF library supports only "
                                                           << (int)NC_MAX_DIMS );
    }
    int nvars;
    fail = nc_inq_nvars( ncFile, &nvars );
    if( nvars > NC_MAX_VARS )
    {
        MB_SET_ERR( MB_FAILURE, "ReadNCDF: File contains " << nvars << " vars but NetCDF library supports only "
                                                           << (int)NC_MAX_VARS );
    }

    // Get the attributes

    // Get the word size, scalar value
    nc_type att_type;
    size_t att_len;
    fail = nc_inq_att( ncFile, NC_GLOBAL, "floating_point_word_size", &att_type, &att_len );
    if( NC_NOERR != fail )
    {
        MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting floating_point_word_size attribute" );
    }
    if( att_type != NC_INT || att_len != 1 )
    {
        MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Word size didn't have type int or size 1" );
    }
    fail = nc_get_att_int( ncFile, NC_GLOBAL, "floating_point_word_size", &IO_WORD_SIZE );
    if( NC_NOERR != fail )
    {
        MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Trouble getting word size" );
    }

    // Exodus version
    fail = nc_inq_att( ncFile, NC_GLOBAL, "version", &att_type, &att_len );
    if( NC_NOERR != fail )
    {
        MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting version attribute" );
    }
    if( att_type != NC_FLOAT || att_len != 1 )
    {
        MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Version didn't have type float or size 1" );
    }
    // float version = temp_att->as_float(0);

    // Read in initial variables
    int temp_dim;
    GET_DIM( temp_dim, "num_dim", numberDimensions_loading );
    GET_DIM( temp_dim, "num_nodes", numberNodes_loading );
    GET_DIM( temp_dim, "num_elem", numberElements_loading );
    GET_DIM( temp_dim, "num_el_blk", numberElementBlocks_loading );
    GET_DIM( temp_dim, "num_fa_blk", numberFaceBlocks_loading );  // for polyhedra blocks
    GET_DIM( temp_dim, "num_elem", numberElements_loading );
    GET_DIM( temp_dim, "num_node_sets", numberNodeSets_loading );
    GET_DIM( temp_dim, "num_side_sets", numberSideSets_loading );
    GET_DIM( temp_dim, "len_string", max_str_length );
    GET_DIM( temp_dim, "len_line", max_line_length );

    // Title; why are we even bothering if we're not going to keep it???
    fail = nc_inq_att( ncFile, NC_GLOBAL, "title", &att_type, &att_len );
    if( NC_NOERR != fail )
    {
        MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting title attribute" );
    }
    if( att_type != NC_CHAR )
    {
        MB_SET_ERR( MB_FAILURE, "ReadNCDF:: title didn't have type char" );
    }
    char* title = new char[att_len + 1];
    fail        = nc_get_att_text( ncFile, NC_GLOBAL, "title", title );
    if( NC_NOERR != fail )
    {
        delete[] title;
        MB_SET_ERR( MB_FAILURE, "ReadNCDF:: trouble getting title" );
    }
    delete[] title;

    return MB_SUCCESS;
}

ErrorCode ReadNCDF::read_nodes( const Tag* file_id_tag )
{
    // Read the nodes into memory

    assert( 0 != ncFile );

    // Create a sequence to hold the node coordinates
    // Get the current number of entities and start at the next slot

    EntityHandle node_handle = 0;
    std::vector< double* > arrays;
    readMeshIface->get_node_coords( 3, numberNodes_loading, MB_START_ID, node_handle, arrays );

    vertexOffset = ID_FROM_HANDLE( node_handle ) - MB_START_ID;

    // Read in the coordinates
    int fail;
    int coord = 0;
    nc_inq_varid( ncFile, "coord", &coord );

    // Single var for all coords
    size_t start[2] = { 0, 0 }, count[2] = { 1, static_cast< size_t >( numberNodes_loading ) };
    if( coord )
    {
        for( int d = 0; d < numberDimensions_loading; ++d )
        {
            start[0] = d;
            fail     = nc_get_vara_double( ncFile, coord, start, count, arrays[d] );
            if( NC_NOERR != fail )
            {
                MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting " << (char)( 'x' + d ) << " coord array" );
            }
        }
    }
    // Var for each coord
    else
    {
        char varname[] = "coord ";
        for( int d = 0; d < numberDimensions_loading; ++d )
        {
            varname[5] = 'x' + (char)d;
            fail       = nc_inq_varid( ncFile, varname, &coord );
            if( NC_NOERR != fail )
            {
                MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting " << (char)( 'x' + d ) << " coord variable" );
            }

            fail = nc_get_vara_double( ncFile, coord, start, &count[1], arrays[d] );
            if( NC_NOERR != fail )
            {
                MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting " << (char)( 'x' + d ) << " coord array" );
            }
        }
    }

    // Zero out any coord values that are in database but not in file
    // (e.g. if MOAB has 3D coords but file is 2D then set Z coords to zero.)
    for( unsigned d = numberDimensions_loading; d < arrays.size(); ++d )
        std::fill( arrays[d], arrays[d] + numberNodes_loading, 0.0 );

    if( file_id_tag )
    {
        Range nodes;
        nodes.insert( node_handle, node_handle + numberNodes_loading - 1 );
        readMeshIface->assign_ids( *file_id_tag, nodes, vertexOffset );
    }

    return MB_SUCCESS;
}

ErrorCode ReadNCDF::read_block_headers( const int* blocks_to_load, const int num_blocks )
{
    // Get the element block ids; keep this in a separate list,
    // which is not offset by blockIdOffset; this list used later for
    // reading block connectivity

    // Get the ids of all the blocks of this file we're reading in
    std::vector< int > block_ids( numberElementBlocks_loading );
    int nc_block_ids = -1;
    GET_1D_INT_VAR( "eb_prop1", nc_block_ids, block_ids );
    if( -1 == nc_block_ids )
    {
        MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting eb_prop1 variable" );
    }

    int exodus_id = 1;

    // If the active_block_id_list is NULL all blocks are active.
    if( NULL == blocks_to_load || 0 == num_blocks )
    {
        blocks_to_load = &block_ids[0];
    }

    std::vector< int > new_blocks( blocks_to_load, blocks_to_load + numberElementBlocks_loading );

    std::vector< int >::iterator iter, end_iter;
    iter     = block_ids.begin();
    end_iter = block_ids.end();

    // Read header information and initialize header-type block information
    int temp_dim;
    std::vector< char > temp_string_storage( max_str_length + 1 );
    char* temp_string = &temp_string_storage[0];
    int block_seq_id  = 1;

    for( ; iter != end_iter; ++iter, block_seq_id++ )
    {
        int num_elements;

        GET_DIMB( temp_dim, temp_string, "num_el_in_blk%d", block_seq_id, num_elements );

        // Don't read element type string for now, since it's an attrib
        // on the connectivity
        // Get the entity type corresponding to this ExoII element type
        // ExoIIElementType elem_type =
        // ExoIIUtil::static_element_name_to_type(element_type_string);

        // Tag each element block(mesh set) with enum for ElementType (SHELL, QUAD4, ....etc)
        ReadBlockData block_data;
        block_data.elemType    = EXOII_MAX_ELEM_TYPE;
        block_data.blockId     = *iter;
        block_data.startExoId  = exodus_id;
        block_data.numElements = num_elements;

        // If block is in 'blocks_to_load'----load it!
        if( std::find( new_blocks.begin(), new_blocks.end(), *iter ) != new_blocks.end() )
            block_data.reading_in = true;
        else
            block_data.reading_in = false;

        blocksLoading.push_back( block_data );
        exodus_id += num_elements;
    }

    return MB_SUCCESS;
}

ErrorCode ReadNCDF::read_face_blocks_headers()
{
    // Get the ids of all the blocks of this file we're reading in
    std::vector< int > fblock_ids( numberFaceBlocks_loading );
    int nc_fblock_ids = -1;
    GET_1D_INT_VAR( "fa_prop1", nc_fblock_ids, fblock_ids );
    if( -1 == nc_fblock_ids )
    {
        MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting fa_prop1 variable" );
    }

    int temp_dim;
    std::vector< char > temp_string_storage( max_str_length + 1 );
    char* temp_string = &temp_string_storage[0];
    // int block_seq_id = 1;
    int exodus_id = 1;

    for( int i = 1; i <= numberFaceBlocks_loading; i++ )
    {
        int num_elements;

        GET_DIMB( temp_dim, temp_string, "num_fa_in_blk%d", i, num_elements );

        // Don't read element type string for now, since it's an attrib
        // on the connectivity
        // Get the entity type corresponding to this ExoII element type
        // ExoIIElementType elem_type =
        // ExoIIUtil::static_element_name_to_type(element_type_string);

        // Tag each element block(mesh set) with enum for ElementType (SHELL, QUAD4, ....etc)
        ReadFaceBlockData fblock_data;
        fblock_data.faceBlockId = i;
        fblock_data.startExoId  = exodus_id;
        fblock_data.numElements = num_elements;

        faceBlocksLoading.push_back( fblock_data );
        exodus_id += num_elements;
    }
    return MB_SUCCESS;
}
ErrorCode ReadNCDF::read_polyhedra_faces()
{
    int temp_dim;
    std::vector< char > temp_string_storage( max_str_length + 1 );
    char* temp_string = &temp_string_storage[0];
    nodesInLoadedBlocks.resize( numberNodes_loading + 1 );
    size_t start[1] = { 0 }, count[1];  // one dim arrays only here!
    std::vector< ReadFaceBlockData >::iterator this_it;
    this_it = faceBlocksLoading.begin();

    int fblock_seq_id = 1;
    std::vector< int > dims;
    int nc_var, fail;

    for( ; this_it != faceBlocksLoading.end(); ++this_it, fblock_seq_id++ )
    {
        int num_fa_in_blk;
        GET_DIMB( temp_dim, temp_string, "num_fa_in_blk%d", fblock_seq_id, num_fa_in_blk );
        int num_nod_per_fa;
        GET_DIMB( temp_dim, temp_string, "num_nod_per_fa%d", fblock_seq_id, num_nod_per_fa );
        // Get the ncdf connect variable and the element type
        INS_ID( temp_string, "fbconn%d", fblock_seq_id );
        GET_VAR( temp_string, nc_var, dims );
        std::vector< int > fbconn;
        fbconn.resize( num_nod_per_fa );
        count[0] = num_nod_per_fa;

        fail = nc_get_vara_int( ncFile, nc_var, start, count, &fbconn[0] );
        if( NC_NOERR != fail )
        {
            MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting face polyhedra connectivity " );
        }
        std::vector< int > fbepecnt;
        fbepecnt.resize( num_fa_in_blk );
        INS_ID( temp_string, "fbepecnt%d", fblock_seq_id );
        GET_VAR( temp_string, nc_var, dims );
        count[0] = num_fa_in_blk;
        fail     = nc_get_vara_int( ncFile, nc_var, start, count, &fbepecnt[0] );
        if( NC_NOERR != fail )
        {
            MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting face polyhedra connectivity " );
        }
        // now we can create some polygons
        std::vector< EntityHandle > polyConn;
        int ix = 0;
        for( int i = 0; i < num_fa_in_blk; i++ )
        {
            polyConn.resize( fbepecnt[i] );
            for( int j = 0; j < fbepecnt[i]; j++ )
            {
                nodesInLoadedBlocks[fbconn[ix]] = 1;
                polyConn[j]                     = vertexOffset + fbconn[ix++];
            }
            EntityHandle newp;
            /*
             *  ErrorCode create_element(const EntityType type,
                                         const EntityHandle *connectivity,
                                         const int num_vertices,
                                         EntityHandle &element_handle)
             */
            if( mdbImpl->create_element( MBPOLYGON, &polyConn[0], fbepecnt[i], newp ) != MB_SUCCESS ) return MB_FAILURE;

            // add the element in order
            polyfaces.push_back( newp );
        }
    }
    return MB_SUCCESS;
}

ErrorCode ReadNCDF::read_elements( const Tag* file_id_tag )
{
    // Read in elements

    int result = 0;

    // Initialize the nodeInLoadedBlocks vector
    if( nodesInLoadedBlocks.size() < (size_t)( numberNodes_loading + 1 ) )
        nodesInLoadedBlocks.resize( numberNodes_loading + 1 );  // this could be repeated?
    memset( &nodesInLoadedBlocks[0], 0, ( numberNodes_loading + 1 ) * sizeof( char ) );

    std::vector< ReadBlockData >::iterator this_it;
    this_it = blocksLoading.begin();

    int temp_dim;
    std::vector< char > temp_string_storage( max_str_length + 1 );
    char* temp_string = &temp_string_storage[0];

    int nc_var;
    int block_seq_id = 1;
    std::vector< int > dims;
    size_t start[2] = { 0, 0 }, count[2];

    for( ; this_it != blocksLoading.end(); ++this_it, block_seq_id++ )
    {
        // If this block isn't to be read in --- continue
        if( !( *this_it ).reading_in ) continue;

        // Get some information about this block
        int block_id       = ( *this_it ).blockId;
        EntityHandle* conn = 0;

        // Get the ncdf connect variable and the element type
        INS_ID( temp_string, "connect%d", block_seq_id );
        GET_VAR( temp_string, nc_var, dims );
        if( -1 == nc_var || 0 == nc_var )
        {  // try other var, for polyhedra blocks
            // it could be polyhedra block, defined by fbconn and NFACED attribute
            INS_ID( temp_string, "facconn%d", block_seq_id );
            GET_VAR( temp_string, nc_var, dims );

            if( -1 == nc_var || 0 == nc_var )
            {
                MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting connec or faccon variable" );
            }
        }
        nc_type att_type;
        size_t att_len;
        int fail = nc_inq_att( ncFile, nc_var, "elem_type", &att_type, &att_len );
        if( NC_NOERR != fail )
        {
            MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting elem type attribute" );
        }
        std::vector< char > dum_str( att_len + 1 );
        fail = nc_get_att_text( ncFile, nc_var, "elem_type", &dum_str[0] );
        if( NC_NOERR != fail )
        {
            MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting elem type" );
        }
        ExoIIElementType elem_type = ExoIIUtil::static_element_name_to_type( &dum_str[0] );
        ( *this_it ).elemType      = elem_type;

        int verts_per_element    = ExoIIUtil::VerticesPerElement[( *this_it ).elemType];
        int number_nodes         = ( *this_it ).numElements * verts_per_element;
        const EntityType mb_type = ExoIIUtil::ExoIIElementMBEntity[( *this_it ).elemType];

        if( mb_type == MBPOLYGON )
        {
            // need to read another variable, num_nod_per_el, and
            int num_nodes_poly_conn;
            GET_DIMB( temp_dim, temp_string, "num_nod_per_el%d", block_seq_id, num_nodes_poly_conn );
            // read the connec1 array, which is of dimensions num_nodes_poly_conn (only one )
            std::vector< int > connec;
            connec.resize( num_nodes_poly_conn );
            count[0] = num_nodes_poly_conn;

            fail = nc_get_vara_int( ncFile, nc_var, start, count, &connec[0] );
            if( NC_NOERR != fail )
            {
                MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting polygon connectivity " );
            }
            // ebepecnt1
            std::vector< int > ebec;
            ebec.resize( this_it->numElements );
            // Get the ncdf connect variable and the element type
            INS_ID( temp_string, "ebepecnt%d", block_seq_id );
            GET_VAR( temp_string, nc_var, dims );
            count[0] = this_it->numElements;
            fail     = nc_get_vara_int( ncFile, nc_var, start, count, &ebec[0] );
            if( NC_NOERR != fail )
            {
                MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting polygon nodes per elem " );
            }
            EntityHandle ms_handle;
            if( mdbImpl->create_meshset( MESHSET_SET | MESHSET_TRACK_OWNER, ms_handle ) != MB_SUCCESS )
                return MB_FAILURE;
            // create polygons one by one, and put them in the list
            // also need to read the number of nodes for each polygon, then create
            std::vector< EntityHandle > polyConn;
            int ix = 0;
            for( int i = 0; i < this_it->numElements; i++ )
            {
                polyConn.resize( ebec[i] );
                for( int j = 0; j < ebec[i]; j++ )
                {
                    nodesInLoadedBlocks[connec[ix]] = 1;
                    polyConn[j]                     = vertexOffset + connec[ix++];
                }
                EntityHandle newp;
                /*
                 *  ErrorCode create_element(const EntityType type,
                                             const EntityHandle *connectivity,
                                             const int num_vertices,
                                             EntityHandle &element_handle)
                 */
                if( mdbImpl->create_element( MBPOLYGON, &polyConn[0], ebec[i], newp ) != MB_SUCCESS ) return MB_FAILURE;

                // add the element in order
                this_it->polys.push_back( newp );
                if( mdbImpl->add_entities( ms_handle, &newp, 1 ) != MB_SUCCESS ) return MB_FAILURE;
            }
            // Set the block id with an offset
            if( mdbImpl->tag_set_data( mMaterialSetTag, &ms_handle, 1, &block_id ) != MB_SUCCESS ) return MB_FAILURE;
            if( mdbImpl->tag_set_data( mGlobalIdTag, &ms_handle, 1, &block_id ) != MB_SUCCESS ) return MB_FAILURE;
        }
        else if( mb_type == MBPOLYHEDRON )
        {
            // need to read another variable, num_fac_per_el
            int num_fac_per_el;
            GET_DIMB( temp_dim, temp_string, "num_fac_per_el%d", block_seq_id, num_fac_per_el );
            // read the fbconn array, which is of dimensions num_nod_per_fa (only one dimension)
            std::vector< int > facconn;
            facconn.resize( num_fac_per_el );
            count[0] = num_fac_per_el;

            fail = nc_get_vara_int( ncFile, nc_var, start, count, &facconn[0] );
            if( NC_NOERR != fail )
            {
                MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting polygon connectivity " );
            }
            // ebepecnt1
            std::vector< int > ebepecnt;
            ebepecnt.resize( this_it->numElements );
            // Get the ncdf connect variable and the element type
            INS_ID( temp_string, "ebepecnt%d", block_seq_id );
            GET_VAR( temp_string, nc_var, dims );
            count[0] = this_it->numElements;
            fail     = nc_get_vara_int( ncFile, nc_var, start, count, &ebepecnt[0] );
            if( NC_NOERR != fail )
            {
                MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting polygon nodes per elem " );
            }
            EntityHandle ms_handle;
            if( mdbImpl->create_meshset( MESHSET_SET | MESHSET_TRACK_OWNER, ms_handle ) != MB_SUCCESS )
                return MB_FAILURE;
            // create polygons one by one, and put them in the list
            // also need to read the number of nodes for each polygon, then create
            std::vector< EntityHandle > polyConn;
            int ix = 0;
            for( int i = 0; i < this_it->numElements; i++ )
            {
                polyConn.resize( ebepecnt[i] );
                for( int j = 0; j < ebepecnt[i]; j++ )
                {
                    polyConn[j] = polyfaces[facconn[ix++] - 1];
                }
                EntityHandle newp;
                /*
                 *  ErrorCode create_element(const EntityType type,
                                             const EntityHandle *connectivity,
                                             const int num_vertices,
                                             EntityHandle &element_handle)
                 */
                if( mdbImpl->create_element( MBPOLYHEDRON, &polyConn[0], ebepecnt[i], newp ) != MB_SUCCESS )
                    return MB_FAILURE;

                // add the element in order
                this_it->polys.push_back( newp );
                if( mdbImpl->add_entities( ms_handle, &newp, 1 ) != MB_SUCCESS ) return MB_FAILURE;
            }
            // Set the block id with an offset
            if( mdbImpl->tag_set_data( mMaterialSetTag, &ms_handle, 1, &block_id ) != MB_SUCCESS ) return MB_FAILURE;
            if( mdbImpl->tag_set_data( mGlobalIdTag, &ms_handle, 1, &block_id ) != MB_SUCCESS ) return MB_FAILURE;
        }
        else  // this is regular
        {
            // Allocate an array to read in connectivity data
            readMeshIface->get_element_connect( this_it->numElements, verts_per_element, mb_type, this_it->startExoId,
                                                this_it->startMBId, conn );

            // Create a range for this sequence of elements
            EntityHandle start_range, end_range;
            start_range = ( *this_it ).startMBId;
            end_range   = start_range + ( *this_it ).numElements - 1;

            Range new_range( start_range, end_range );
            // Range<EntityHandle> new_range((*this_it).startMBId,
            //                              (*this_it).startMBId + (*this_it).numElements - 1);

            // Create a MBSet for this block and set the material tag

            EntityHandle ms_handle;
            if( mdbImpl->create_meshset( MESHSET_SET | MESHSET_TRACK_OWNER, ms_handle ) != MB_SUCCESS )
                return MB_FAILURE;

            if( mdbImpl->add_entities( ms_handle, new_range ) != MB_SUCCESS ) return MB_FAILURE;

            int mid_nodes[4];
            CN::HasMidNodes( mb_type, verts_per_element, mid_nodes );
            if( mdbImpl->tag_set_data( mHasMidNodesTag, &ms_handle, 1, mid_nodes ) != MB_SUCCESS ) return MB_FAILURE;

            // Just a check because the following code won't work if this case fails
            assert( sizeof( EntityHandle ) >= sizeof( int ) );

            // tmp_ptr is of type int* and points at the same place as conn
            int* tmp_ptr = reinterpret_cast< int* >( conn );

            // Read the connectivity into that memory,  this will take only part of the array
            // 1/2 if sizeof(EntityHandle) == 64 bits.
            count[0] = this_it->numElements;
            count[1] = verts_per_element;
            fail     = nc_get_vara_int( ncFile, nc_var, start, count, tmp_ptr );
            if( NC_NOERR != fail )
            {
                MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting connectivity" );
            }

            // Convert from exodus indices to vertex handles.
            // Iterate backwards in case handles are larger than ints.
            for( int i = number_nodes - 1; i >= 0; --i )
            {
                if( (unsigned)tmp_ptr[i] >= nodesInLoadedBlocks.size() )
                {
                    MB_SET_ERR( MB_FAILURE, "Invalid node ID in block connectivity" );
                }
                nodesInLoadedBlocks[tmp_ptr[i]] = 1;
                conn[i]                         = static_cast< EntityHandle >( tmp_ptr[i] ) + vertexOffset;
            }

            // Adjust connectivity order if necessary
            const int* reorder = exodus_elem_order_map[mb_type][verts_per_element];
            if( reorder ) ReadUtilIface::reorder( reorder, conn, this_it->numElements, verts_per_element );

            readMeshIface->update_adjacencies( ( *this_it ).startMBId, ( *this_it ).numElements,
                                               ExoIIUtil::VerticesPerElement[( *this_it ).elemType], conn );

            if( result == -1 )
            {
                MB_SET_ERR( MB_FAILURE, "ReadNCDF:: error getting element connectivity for block " << block_id );
            }

            // Set the block id with an offset
            if( mdbImpl->tag_set_data( mMaterialSetTag, &ms_handle, 1, &block_id ) != MB_SUCCESS ) return MB_FAILURE;
            if( mdbImpl->tag_set_data( mGlobalIdTag, &ms_handle, 1, &block_id ) != MB_SUCCESS ) return MB_FAILURE;

            if( file_id_tag )
            {
                Range range;
                range.insert( this_it->startMBId, this_it->startMBId + this_it->numElements - 1 );
                readMeshIface->assign_ids( *file_id_tag, range, this_it->startExoId );
            }
        }  // end regular block (not polygon)
    }

    return MB_SUCCESS;
}

ErrorCode ReadNCDF::read_global_ids()
{
    // Read in the map from the exodus file
    std::vector< int > ids( std::max( numberElements_loading, numberNodes_loading ) );

    int varid = -1;
    GET_1D_INT_VAR( "elem_map", varid, ids );
    if( -1 != varid )
    {
        std::vector< ReadBlockData >::iterator iter;
        int id_pos = 0;
        for( iter = blocksLoading.begin(); iter != blocksLoading.end(); ++iter )
        {
            if( iter->reading_in )
            {
                if( iter->elemType != EXOII_POLYGON && iter->elemType != EXOII_POLYHEDRON )
                {
                    if( iter->startMBId != 0 )
                    {
                        Range range( iter->startMBId, iter->startMBId + iter->numElements - 1 );
                        ErrorCode error = mdbImpl->tag_set_data( mGlobalIdTag, range, &ids[id_pos] );
                        if( error != MB_SUCCESS ) return error;
                        id_pos += iter->numElements;
                    }
                    else
                        return MB_FAILURE;
                }
                else  // polygons or polyhedrons; use id from elements
                {
                    for( std::vector< EntityHandle >::iterator eit = iter->polys.begin(); eit != iter->polys.end();
                         eit++, id_pos++ )
                    {
                        EntityHandle peh = *eit;
                        if( mdbImpl->tag_set_data( mGlobalIdTag, &peh, 1, &ids[id_pos] ) != MB_SUCCESS )
                            return MB_FAILURE;
                    }
                }
            }
        }
    }

    // Read in node map next
    varid = -1;
    GET_1D_INT_VAR( "node_num_map", varid, ids );
    if( -1 != varid )
    {
        Range range( MB_START_ID + vertexOffset, MB_START_ID + vertexOffset + numberNodes_loading - 1 );
        ErrorCode error = mdbImpl->tag_set_data( mGlobalIdTag, range, &ids[0] );
        if( MB_SUCCESS != error ) MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem setting node global ids" );
    }

    return MB_SUCCESS;
}

ErrorCode ReadNCDF::read_nodesets()
{
    // Read in the nodesets for the model

    if( 0 == numberNodeSets_loading ) return MB_SUCCESS;
    std::vector< int > id_array( numberNodeSets_loading );

    // Read in the nodeset ids
    int nc_var;
    GET_1D_INT_VAR( "ns_prop1", nc_var, id_array );
    if( -1 == nc_var )
    {
        MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting ns_prop1 variable" );
    }

    // Use a vector of ints to read node handles
    std::vector< int > node_handles;

    int i;
    std::vector< char > temp_string_storage( max_str_length + 1 );
    char* temp_string = &temp_string_storage[0];
    int temp_dim;
    for( i = 0; i < numberNodeSets_loading; i++ )
    {
        // Get nodeset parameters
        int number_nodes_in_set;
        int number_dist_factors_in_set;

        GET_DIMB( temp_dim, temp_string, "num_nod_ns%d", i + 1, number_nodes_in_set );
        GET_DIMB( temp_dim, temp_string, "num_df_ns%d", i + 1, number_dist_factors_in_set );

        // Need to new a vector to store dist. factors
        // This vector * gets stored as a tag on the sideset meshset
        std::vector< double > temp_dist_factor_vector( number_nodes_in_set );
        if( number_dist_factors_in_set != 0 )
        {
            INS_ID( temp_string, "dist_fact_ns%d", i + 1 );
            GET_1D_DBL_VAR( temp_string, temp_dim, temp_dist_factor_vector );
            if( -1 == temp_dim )
            {
                MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting dist fact variable" );
            }
        }

        // Size new arrays and get ids and distribution factors
        if( node_handles.size() < (unsigned int)number_nodes_in_set )
        {
            node_handles.reserve( number_nodes_in_set );
            node_handles.resize( number_nodes_in_set );
        }

        INS_ID( temp_string, "node_ns%d", i + 1 );
        int temp_var = -1;
        GET_1D_INT_VAR( temp_string, temp_var, node_handles );
        if( -1 == temp_var )
        {
            MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting nodeset node variable" );
        }

        // Maybe there is already a nodesets meshset here we can append to
        Range child_meshsets;
        if( mdbImpl->get_entities_by_handle( 0, child_meshsets ) != MB_SUCCESS ) return MB_FAILURE;

        child_meshsets = subtract( child_meshsets, initRange );

        Range::iterator iter, end_iter;
        iter     = child_meshsets.begin();
        end_iter = child_meshsets.end();

        EntityHandle ns_handle = 0;
        for( ; iter != end_iter; ++iter )
        {
            int nodeset_id;
            if( mdbImpl->tag_get_data( mDirichletSetTag, &( *iter ), 1, &nodeset_id ) != MB_SUCCESS ) continue;

            if( id_array[i] == nodeset_id )
            {
                // Found the meshset
                ns_handle = *iter;
                break;
            }
        }

        std::vector< EntityHandle > nodes_of_nodeset;
        if( ns_handle )
            if( mdbImpl->get_entities_by_handle( ns_handle, nodes_of_nodeset, true ) != MB_SUCCESS ) return MB_FAILURE;

        // Make these into entity handles
        // TODO: could we have read it into EntityHandle sized array in the first place?
        int j, temp;
        std::vector< EntityHandle > nodes;
        std::vector< double > dist_factor_vector;
        for( j = 0; j < number_nodes_in_set; j++ )
        {
            // See if this node is one we're currently reading in
            if( nodesInLoadedBlocks[node_handles[j]] == 1 )
            {
                // Make sure that it already isn't in a nodeset
                unsigned int node_id = CREATE_HANDLE( MBVERTEX, node_handles[j] + vertexOffset, temp );
                if( !ns_handle ||
                    std::find( nodes_of_nodeset.begin(), nodes_of_nodeset.end(), node_id ) == nodes_of_nodeset.end() )
                {
                    nodes.push_back( node_id );

                    if( number_dist_factors_in_set != 0 ) dist_factor_vector.push_back( temp_dist_factor_vector[j] );
                }
            }
        }

        // No nodes to add
        if( nodes.empty() ) continue;

        // If there was no meshset found --> create one
        if( ns_handle == 0 )
        {
            if( mdbImpl->create_meshset( MESHSET_ORDERED | MESHSET_TRACK_OWNER, ns_handle ) != MB_SUCCESS )
                return MB_FAILURE;

            // Set a tag signifying dirichlet bc
            // TODO: create this tag another way

            int nodeset_id = id_array[i];
            if( mdbImpl->tag_set_data( mDirichletSetTag, &ns_handle, 1, &nodeset_id ) != MB_SUCCESS ) return MB_FAILURE;
            if( mdbImpl->tag_set_data( mGlobalIdTag, &ns_handle, 1, &nodeset_id ) != MB_SUCCESS ) return MB_FAILURE;

            if( !dist_factor_vector.empty() )
            {
                int size         = dist_factor_vector.size();
                const void* data = &dist_factor_vector[0];
                if( mdbImpl->tag_set_by_ptr( mDistFactorTag, &ns_handle, 1, &data, &size ) != MB_SUCCESS )
                    return MB_FAILURE;
            }
        }
        else if( !dist_factor_vector.empty() )
        {
            // Append dist factors to vector
            const void* ptr = 0;
            int size        = 0;
            if( mdbImpl->tag_get_by_ptr( mDistFactorTag, &ns_handle, 1, &ptr, &size ) != MB_SUCCESS ) return MB_FAILURE;

            const double* data = reinterpret_cast< const double* >( ptr );
            dist_factor_vector.reserve( dist_factor_vector.size() + size );
            std::copy( data, data + size, std::back_inserter( dist_factor_vector ) );
            size = dist_factor_vector.size();
            ptr  = &dist_factor_vector[0];
            if( mdbImpl->tag_set_by_ptr( mDistFactorTag, &ns_handle, 1, &ptr, &size ) != MB_SUCCESS ) return MB_FAILURE;
        }

        // Add the nodes to the meshset
        if( mdbImpl->add_entities( ns_handle, &nodes[0], nodes.size() ) != MB_SUCCESS ) return MB_FAILURE;
    }

    return MB_SUCCESS;
}

ErrorCode ReadNCDF::read_sidesets()
{
    // Uhhh if you read the same file (previously_read==true) then blow away the
    // sidesets pertaining to this file?  How do you do that? If you read a
    // new file make sure all the offsets are correct.

    // If not loading any sidesets -- exit
    if( 0 == numberSideSets_loading ) return MB_SUCCESS;

    // Read in the sideset ids
    std::vector< int > id_array( numberSideSets_loading );
    int temp_var = -1;
    GET_1D_INT_VAR( "ss_prop1", temp_var, id_array );
    if( -1 == temp_var )
    {
        MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting ss_prop1 variable" );
    }

    // Create a side set for each one
    int number_sides_in_set;
    int number_dist_factors_in_set;

    // Maybe there is already a sidesets meshset here we can append to
    Range child_meshsets;
    if( mdbImpl->get_entities_by_type( 0, MBENTITYSET, child_meshsets ) != MB_SUCCESS ) return MB_FAILURE;

    child_meshsets = subtract( child_meshsets, initRange );

    Range::iterator iter, end_iter;

    int i;
    std::vector< char > temp_string_storage( max_str_length + 1 );
    char* temp_string = &temp_string_storage[0];
    int temp_dim;
    for( i = 0; i < numberSideSets_loading; i++ )
    {
        // Get sideset parameters
        GET_DIMB( temp_dim, temp_string, "num_side_ss%d", i + 1, number_sides_in_set );
        GET_DIMB( temp_dim, temp_string, "num_df_ss%d", i + 1, number_dist_factors_in_set );

        // Size new arrays and get element and side lists
        std::vector< int > side_list( number_sides_in_set );
        std::vector< int > element_list( number_sides_in_set );
        INS_ID( temp_string, "side_ss%d", i + 1 );
        temp_var = -1;
        GET_1D_INT_VAR( temp_string, temp_var, side_list );
        if( -1 == temp_var )
        {
            MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting sideset side variable" );
        }

        INS_ID( temp_string, "elem_ss%d", i + 1 );
        temp_var = -1;
        GET_1D_INT_VAR( temp_string, temp_var, element_list );
        if( -1 == temp_var )
        {
            MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting sideset elem variable" );
        }

        std::vector< double > temp_dist_factor_vector;
        std::vector< EntityHandle > entities_to_add, reverse_entities;
        // Create the sideset entities
        if( create_ss_elements( &element_list[0], &side_list[0], number_sides_in_set, number_dist_factors_in_set,
                                entities_to_add, reverse_entities, temp_dist_factor_vector, i + 1 ) != MB_SUCCESS )
            return MB_FAILURE;

        // If there are elements to add
        if( !entities_to_add.empty() || !reverse_entities.empty() )
        {
            iter     = child_meshsets.begin();
            end_iter = child_meshsets.end();

            EntityHandle ss_handle = 0;
            for( ; iter != end_iter; ++iter )
            {
                int sideset_id;
                if( mdbImpl->tag_get_data( mNeumannSetTag, &( *iter ), 1, &sideset_id ) != MB_SUCCESS ) continue;

                if( id_array[i] == sideset_id )
                {
                    // Found the meshset
                    ss_handle = *iter;
                    break;
                }
            }

            // If we didn't find a sideset already
            if( ss_handle == 0 )
            {
                if( mdbImpl->create_meshset( MESHSET_ORDERED | MESHSET_TRACK_OWNER, ss_handle ) != MB_SUCCESS )
                    return MB_FAILURE;

                if( ss_handle == 0 ) return MB_FAILURE;

                int sideset_id = id_array[i];
                if( mdbImpl->tag_set_data( mNeumannSetTag, &ss_handle, 1, &sideset_id ) != MB_SUCCESS )
                    return MB_FAILURE;
                if( mdbImpl->tag_set_data( mGlobalIdTag, &ss_handle, 1, &sideset_id ) != MB_SUCCESS ) return MB_FAILURE;

                if( !reverse_entities.empty() )
                {
                    // Also make a reverse set to put in this set
                    EntityHandle reverse_set;
                    if( mdbImpl->create_meshset( MESHSET_SET | MESHSET_TRACK_OWNER, reverse_set ) != MB_SUCCESS )
                        return MB_FAILURE;

                    // Add the reverse set to the sideset set and the entities to the reverse set
                    ErrorCode result = mdbImpl->add_entities( ss_handle, &reverse_set, 1 );
                    if( MB_SUCCESS != result ) return result;

                    result = mdbImpl->add_entities( reverse_set, &reverse_entities[0], reverse_entities.size() );
                    if( MB_SUCCESS != result ) return result;

                    // Set the reverse tag
                    Tag sense_tag;
                    int dum_sense = 0;
                    result        = mdbImpl->tag_get_handle( "NEUSET_SENSE", 1, MB_TYPE_INTEGER, sense_tag,
                                                             MB_TAG_SPARSE | MB_TAG_CREAT, &dum_sense );
                    if( result != MB_SUCCESS ) return result;
                    dum_sense = -1;
                    result    = mdbImpl->tag_set_data( sense_tag, &reverse_set, 1, &dum_sense );
                    if( result != MB_SUCCESS ) return result;
                }
            }

            if( mdbImpl->add_entities( ss_handle, ( entities_to_add.empty() ) ? NULL : &entities_to_add[0],
                                       entities_to_add.size() ) != MB_SUCCESS )
                return MB_FAILURE;

            // Distribution factor stuff
            if( number_dist_factors_in_set )
            {
                // If this sideset does not already has a distribution factor array...set one
                const void* ptr = 0;
                int size        = 0;
                if( MB_SUCCESS == mdbImpl->tag_get_by_ptr( mDistFactorTag, &ss_handle, 1, &ptr, &size ) )
                {
                    const double* data = reinterpret_cast< const double* >( ptr );
                    std::copy( data, data + size, std::back_inserter( temp_dist_factor_vector ) );
                }

                ptr  = &temp_dist_factor_vector[0];
                size = temp_dist_factor_vector.size();
                if( mdbImpl->tag_set_by_ptr( mDistFactorTag, &ss_handle, 1, &ptr, &size ) != MB_SUCCESS )
                    return MB_FAILURE;
            }
        }
    }

    return MB_SUCCESS;
}

ErrorCode ReadNCDF::create_ss_elements( int* element_ids,
                                        int* side_list,
                                        int num_sides,
                                        int num_dist_factors,
                                        std::vector< EntityHandle >& entities_to_add,
                                        std::vector< EntityHandle >& reverse_entities,
                                        std::vector< double >& dist_factor_vector,
                                        int ss_seq_id )
{
    // Determine entity type from element id

    // If there are dist. factors, create a vector to hold the array
    // and place this array as a tag onto the sideset meshset

    std::vector< double > temp_dist_factor_vector( num_dist_factors );
    std::vector< char > temp_string_storage( max_str_length + 1 );
    char* temp_string = &temp_string_storage[0];
    int temp_var;
    if( num_dist_factors )
    {
        INS_ID( temp_string, "dist_fact_ss%d", ss_seq_id );
        temp_var = -1;
        GET_1D_DBL_VAR( temp_string, temp_var, temp_dist_factor_vector );
        if( -1 == temp_var )
        {
            MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting dist fact variable" );
        }
    }

    EntityType subtype;
    int num_side_nodes, num_elem_nodes;
    const EntityHandle* nodes;
    std::vector< EntityHandle > connectivity;
    int side_node_idx[32];

    int df_index = 0;
    int sense    = 0;
    for( int i = 0; i < num_sides; i++ )
    {
        ExoIIElementType exoii_type;
        ReadBlockData block_data;
        block_data.elemType = EXOII_MAX_ELEM_TYPE;

        if( find_side_element_type( element_ids[i], exoii_type, block_data, df_index, side_list[i] ) != MB_SUCCESS )
            continue;  // Isn't being read in this time

        EntityType type = ExoIIUtil::ExoIIElementMBEntity[exoii_type];

        EntityHandle ent_handle = element_ids[i] - block_data.startExoId + block_data.startMBId;

        const int side_num = side_list[i] - 1;
        if( type == MBHEX )
        {
            // Get the nodes of the element
            if( mdbImpl->get_connectivity( ent_handle, nodes, num_elem_nodes ) != MB_SUCCESS ) return MB_FAILURE;

            CN::SubEntityNodeIndices( type, num_elem_nodes, 2, side_num, subtype, num_side_nodes, side_node_idx );
            if( num_side_nodes <= 0 ) return MB_FAILURE;

            connectivity.resize( num_side_nodes );
            for( int k = 0; k < num_side_nodes; ++k )
                connectivity[k] = nodes[side_node_idx[k]];

            if( MB_SUCCESS != create_sideset_element( connectivity, subtype, ent_handle, sense ) ) return MB_FAILURE;
            if( 1 == sense )
                entities_to_add.push_back( ent_handle );
            else
                reverse_entities.push_back( ent_handle );

            // Read in distribution factor array
            if( num_dist_factors )
            {
                for( int k = 0; k < 4; k++ )
                    dist_factor_vector.push_back( temp_dist_factor_vector[df_index++] );
            }
        }

        // If it is a Tet
        else if( type == MBTET )
        {
            // Get the nodes of the element
            if( mdbImpl->get_connectivity( ent_handle, nodes, num_elem_nodes ) != MB_SUCCESS ) return MB_FAILURE;

            CN::SubEntityNodeIndices( type, num_elem_nodes, 2, side_num, subtype, num_side_nodes, side_node_idx );
            if( num_side_nodes <= 0 ) return MB_FAILURE;

            connectivity.resize( num_side_nodes );
            for( int k = 0; k < num_side_nodes; ++k )
                connectivity[k] = nodes[side_node_idx[k]];

            if( MB_SUCCESS != create_sideset_element( connectivity, subtype, ent_handle, sense ) ) return MB_FAILURE;
            if( 1 == sense )
                entities_to_add.push_back( ent_handle );
            else
                reverse_entities.push_back( ent_handle );

            // Read in distribution factor array
            if( num_dist_factors )
            {
                for( int k = 0; k < 3; k++ )
                    dist_factor_vector.push_back( temp_dist_factor_vector[df_index++] );
            }
        }
        else if( type == MBQUAD && exoii_type >= EXOII_SHELL && exoii_type <= EXOII_SHELL9 )
        {
            // ent_handle = CREATE_HANDLE(MBQUAD, base_id, error);

            // Just use this quad
            if( side_list[i] == 1 )
            {
                if( 1 == sense )
                    entities_to_add.push_back( ent_handle );
                else
                    reverse_entities.push_back( ent_handle );

                if( num_dist_factors )
                {
                    for( int k = 0; k < 4; k++ )
                        dist_factor_vector.push_back( temp_dist_factor_vector[df_index++] );
                }

                continue;
            }
            else if( side_list[i] == 2 )
            {
                reverse_entities.push_back( ent_handle );

                if( num_dist_factors )
                {
                    for( int k = 0; k < 4; k++ )
                        dist_factor_vector.push_back( temp_dist_factor_vector[df_index++] );
                }

                continue;
            }
            else
            {
                // Get the nodes of the element
                if( mdbImpl->get_connectivity( ent_handle, nodes, num_elem_nodes ) != MB_SUCCESS ) return MB_FAILURE;

                CN::SubEntityNodeIndices( type, num_elem_nodes, 1, side_num - 2, subtype, num_side_nodes,
                                          side_node_idx );
                if( num_side_nodes <= 0 ) return MB_FAILURE;

                connectivity.resize( num_side_nodes );
                for( int k = 0; k < num_side_nodes; ++k )
                    connectivity[k] = nodes[side_node_idx[k]];

                if( MB_SUCCESS != create_sideset_element( connectivity, subtype, ent_handle, sense ) )
                    return MB_FAILURE;
                if( 1 == sense )
                    entities_to_add.push_back( ent_handle );
                else
                    reverse_entities.push_back( ent_handle );

                if( num_dist_factors )
                {
                    for( int k = 0; k < 2; k++ )
                        dist_factor_vector.push_back( temp_dist_factor_vector[df_index++] );
                }
            }
        }
        // If it is a Quad
        else if( type == MBQUAD )
        {
            // Get the nodes of the element
            if( mdbImpl->get_connectivity( ent_handle, nodes, num_elem_nodes ) != MB_SUCCESS ) return MB_FAILURE;

            CN::SubEntityNodeIndices( type, num_elem_nodes, 1, side_num, subtype, num_side_nodes, side_node_idx );
            if( num_side_nodes <= 0 ) return MB_FAILURE;

            connectivity.resize( num_side_nodes );
            for( int k = 0; k < num_side_nodes; ++k )
                connectivity[k] = nodes[side_node_idx[k]];

            if( MB_SUCCESS != create_sideset_element( connectivity, subtype, ent_handle, sense ) ) return MB_FAILURE;
            if( 1 == sense )
                entities_to_add.push_back( ent_handle );
            else
                reverse_entities.push_back( ent_handle );

            // Read in distribution factor array
            if( num_dist_factors )
            {
                for( int k = 0; k < 2; k++ )
                    dist_factor_vector.push_back( temp_dist_factor_vector[df_index++] );
            }
        }
        else if( type == MBTRI )
        {
            int side_offset = 0;
            if( number_dimensions() == 3 && side_list[i] <= 2 )
            {
                if( 1 == sense )
                    entities_to_add.push_back( ent_handle );
                else
                    reverse_entities.push_back( ent_handle );
                if( num_dist_factors )
                {
                    for( int k = 0; k < 3; k++ )
                        dist_factor_vector.push_back( temp_dist_factor_vector[df_index++] );
                }
            }
            else
            {
                if( number_dimensions() == 3 )
                {
                    if( side_list[i] > 2 ) side_offset = 2;
                }

                // Get the nodes of the element
                if( mdbImpl->get_connectivity( ent_handle, nodes, num_elem_nodes ) != MB_SUCCESS ) return MB_FAILURE;

                CN::SubEntityNodeIndices( type, num_elem_nodes, 1, side_num - side_offset, subtype, num_side_nodes,
                                          side_node_idx );
                if( num_side_nodes <= 0 ) return MB_FAILURE;

                connectivity.resize( num_side_nodes );
                for( int k = 0; k < num_side_nodes; ++k )
                    connectivity[k] = nodes[side_node_idx[k]];

                if( MB_SUCCESS != create_sideset_element( connectivity, subtype, ent_handle, sense ) )
                    return MB_FAILURE;
                if( 1 == sense )
                    entities_to_add.push_back( ent_handle );
                else
                    reverse_entities.push_back( ent_handle );

                if( num_dist_factors )
                {
                    for( int k = 0; k < 2; k++ )
                        dist_factor_vector.push_back( temp_dist_factor_vector[df_index++] );
                }
            }
        }
    }

    return MB_SUCCESS;
}

ErrorCode ReadNCDF::create_sideset_element( const std::vector< EntityHandle >& connectivity,
                                            EntityType type,
                                            EntityHandle& handle,
                                            int& sense )
{
    // Get adjacent entities
    ErrorCode error = MB_SUCCESS;
    int to_dim      = CN::Dimension( type );
    // for matching purposes , consider only corner nodes
    // basically, assume everything is matching if the corner nodes are matching, and
    // the number of total nodes is the same
    int nb_corner_nodes = CN::VerticesPerEntity( type );
    std::vector< EntityHandle > adj_ent;
    mdbImpl->get_adjacencies( &( connectivity[0] ), 1, to_dim, false, adj_ent );

    // For each entity, see if we can find a match
    // If we find a match, return it
    bool match_found = false;
    std::vector< EntityHandle > match_conn;
    // By default, assume positive sense
    sense = 1;
    for( unsigned int i = 0; i < adj_ent.size() && match_found == false; i++ )
    {
        // Get the connectivity
        error = mdbImpl->get_connectivity( &( adj_ent[i] ), 1, match_conn );
        if( error != MB_SUCCESS ) continue;

        // Make sure they have the same number of vertices (higher order elements ?)
        if( match_conn.size() != connectivity.size() ) continue;

        // Find a matching node
        std::vector< EntityHandle >::iterator iter;
        std::vector< EntityHandle >::iterator end_corner = match_conn.begin() + nb_corner_nodes;
        iter                                             = std::find( match_conn.begin(), end_corner, connectivity[0] );
        if( iter == end_corner ) continue;

        // Rotate to match connectivity
        // rotate only corner nodes, ignore the rest from now on
        std::rotate( match_conn.begin(), iter, end_corner );

        bool they_match = true;
        for( int j = 1; j < nb_corner_nodes; j++ )
        {
            if( connectivity[j] != match_conn[j] )
            {
                they_match = false;
                break;
            }
        }

        // If we didn't get a match
        if( !they_match )
        {
            // Try the opposite sense
            they_match = true;

            int k = nb_corner_nodes - 1;
            for( int j = 1; j < nb_corner_nodes; )
            {
                if( connectivity[j] != match_conn[k] )
                {
                    they_match = false;
                    break;
                }
                ++j;
                --k;
            }
            // If they matched here, sense is reversed
            if( they_match ) sense = -1;
        }
        match_found = they_match;
        if( match_found ) handle = adj_ent[i];
    }

    // If we didn't find a match, create an element
    if( !match_found ) error = mdbImpl->create_element( type, &connectivity[0], connectivity.size(), handle );

    return error;
}

ErrorCode ReadNCDF::find_side_element_type( const int exodus_id,
                                            ExoIIElementType& elem_type,
                                            ReadBlockData& block_data,
                                            int& df_index,
                                            int side_id )
{
    std::vector< ReadBlockData >::iterator iter, end_iter;
    iter      = blocksLoading.begin();
    end_iter  = blocksLoading.end();
    elem_type = EXOII_MAX_ELEM_TYPE;

    for( ; iter != end_iter; ++iter )
    {
        if( exodus_id >= ( *iter ).startExoId && exodus_id < ( *iter ).startExoId + ( *iter ).numElements )
        {
            elem_type = ( *iter ).elemType;

            // If we're not reading this block in
            if( !( *iter ).reading_in )
            {
                // Offset df_indes according to type
                if( elem_type >= EXOII_HEX && elem_type <= EXOII_HEX27 )
                    df_index += 4;
                else if( elem_type >= EXOII_TETRA && elem_type <= EXOII_TETRA14 )
                    df_index += 3;
                else if( elem_type >= EXOII_QUAD && elem_type <= EXOII_QUAD9 )
                    df_index += 2;
                else if( elem_type >= EXOII_SHELL && elem_type <= EXOII_SHELL9 )
                {
                    if( side_id == 1 || side_id == 2 )
                        df_index += 4;
                    else
                        df_index += 2;
                }
                else if( elem_type >= EXOII_TRI && elem_type <= EXOII_TRI7 )
                    df_index += 3;

                return MB_FAILURE;
            }

            block_data = *iter;

            return MB_SUCCESS;
        }
    }
    return MB_FAILURE;
}

ErrorCode ReadNCDF::read_qa_records( EntityHandle file_set )
{
    std::vector< std::string > qa_records;
    read_qa_information( qa_records );

    std::vector< char > tag_data;
    for( std::vector< std::string >::iterator i = qa_records.begin(); i != qa_records.end(); ++i )
    {
        std::copy( i->begin(), i->end(), std::back_inserter( tag_data ) );
        tag_data.push_back( '\0' );
    }

    // If there were qa_records...tag them to the mCurrentMeshHandle
    if( !tag_data.empty() )
    {
        const void* ptr = &tag_data[0];
        int size        = tag_data.size();
        if( mdbImpl->tag_set_by_ptr( mQaRecordTag, &file_set, 1, &ptr, &size ) != MB_SUCCESS )
        {
            return MB_FAILURE;
        }
    }

    return MB_SUCCESS;
}

ErrorCode ReadNCDF::read_qa_information( std::vector< std::string >& qa_record_list )
{
    // Inquire on the genesis file to find the number of qa records

    int number_records = 0;

    int temp_dim;
    GET_DIM( temp_dim, "num_qa_rec", number_records );
    std::vector< char > data( max_str_length + 1 );

    for( int i = 0; i < number_records; i++ )
    {
        for( int j = 0; j < 4; j++ )
        {
            data[max_str_length] = '\0';
            if( read_qa_string( &data[0], i, j ) != MB_SUCCESS ) return MB_FAILURE;

            qa_record_list.push_back( &data[0] );
        }
    }

    return MB_SUCCESS;
}

ErrorCode ReadNCDF::read_qa_string( char* temp_string, int record_number, int record_position )
{
    std::vector< int > dims;
    int temp_var = -1;
    GET_VAR( "qa_records", temp_var, dims );
    if( -1 == temp_var )
    {
        MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting qa record variable" );
        return MB_FAILURE;
    }
    size_t count[3], start[3];
    start[0] = record_number;
    start[1] = record_position;
    start[2] = 0;
    count[0] = 1;
    count[1] = 1;
    count[2] = max_str_length;
    int fail = nc_get_vara_text( ncFile, temp_var, start, count, temp_string );
    if( NC_NOERR != fail )
    {
        MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem setting current record number variable position" );
    }
    // Get the variable id in the exodus file

    return MB_SUCCESS;
}

// The cub_file_set contains the mesh to be updated. There could exist other
// file sets that should be kept separate, such as the geometry file set from
// ReadCGM.
ErrorCode ReadNCDF::update( const char* exodus_file_name,
                            const FileOptions& opts,
                            const int num_blocks,
                            const int* blocks_to_load,
                            const EntityHandle cub_file_set )
{
    // Function : updating current database from new exodus_file.
    // Creator:   Jane Hu
    // opts is currently designed as following
    // tdata = <var_name>[, time][,op][,destination]
    // where var_name show the tag name to be updated, this version just takes
    // coord.
    // time is the optional, and it gives time step of each of the mesh
    // info in exodus file. It start from 1.
    // op is the operation that is going to be performed on the var_name info.
    // currently support 'set'
    // destination shows where to store the updated info, currently assume it is
    // stored in the same database by replacing the old info if there's no input
    // for destination, or the destination data is given in exodus format and just
    // need to update the coordinates.
    // Assumptions:
    // 1. Assume the num_el_blk's in both DB and update exodus file are the same.
    // 2. Assume num_el_in_blk1...num_el_in_blk(num_el_blk) numbers are matching, may in
    // different order. example: num_el_in_blk11 = num_el_in_blk22 && num_el_in_blk12 =
    // num_el_in_blk21.
    // 3. In exodus file, get node_num_map
    // 4. loop through the node_num_map, use it to find the node in the cub file.
    // 5. Replace coord[0][n] with coordx[m] + vals_nod_var1(time_step, m) for all directions for
    // matching nodes.
    //  Test: try to get hexes

    // *******************************************************************
    // Move nodes to their deformed locations.
    // *******************************************************************
    ErrorCode rval;
    std::string s;
    rval = opts.get_str_option( "tdata", s );
    if( MB_SUCCESS != rval )
    {
        MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem reading file options" );
    }
    std::vector< std::string > tokens;
    tokenize( s, tokens, "," );

    // 1. Check for time step to find the match time
    int time_step = 1;
    if( tokens.size() > 1 && !tokens[1].empty() )
    {
        const char* time_s = tokens[1].c_str();
        char* endptr;
        long int pval   = strtol( time_s, &endptr, 0 );
        std::string end = endptr;
        if( !end.empty() )  // Syntax error
            return MB_TYPE_OUT_OF_RANGE;

        // Check for overflow (parsing long int, returning int)
        time_step = pval;
        if( pval != (long int)time_step ) return MB_TYPE_OUT_OF_RANGE;
        if( time_step <= 0 ) return MB_TYPE_OUT_OF_RANGE;
    }

    // 2. Check for the operations, currently support set.
    const char* op;
    if( tokens.size() < 3 || ( tokens[2] != "set" && tokens[2] != "add" ) )
    {
        MB_SET_ERR( MB_TYPE_OUT_OF_RANGE, "ReadNCDF: invalid operation specified for update" );
    }
    op = tokens[2].c_str();

    // Open netcdf/exodus file
    ncFile   = 0;
    int fail = nc_open( exodus_file_name, 0, &ncFile );
    if( !ncFile )
    {
        MB_SET_ERR( MB_FILE_DOES_NOT_EXIST, "ReadNCDF:: problem opening Netcdf/Exodus II file " << exodus_file_name );
    }

    rval = read_exodus_header();MB_CHK_ERR( rval );

    // Check to make sure that the requested time step exists
    int ncdim = -1;
    int max_time_steps;
    GET_DIM( ncdim, "time_step", max_time_steps );
    if( -1 == ncdim )
    {
        std::cout << "ReadNCDF: could not get number of time steps" << std::endl;
        return MB_FAILURE;
    }
    std::cout << "  Maximum time step=" << max_time_steps << std::endl;
    if( max_time_steps < time_step )
    {
        std::cout << "ReadNCDF: time step is greater than max_time_steps" << std::endl;
        return MB_FAILURE;
    }

    // Get the time
    std::vector< double > times( max_time_steps );
    int nc_var = -1;
    GET_1D_DBL_VAR( "time_whole", nc_var, times );
    if( -1 == nc_var )
    {
        std::cout << "ReadNCDF: unable to get time variable" << std::endl;
    }
    else
    {
        std::cout << "  Step " << time_step << " is at " << times[time_step - 1] << " seconds" << std::endl;
    }

    // Read in the node_num_map.
    std::vector< int > ptr( numberNodes_loading );

    int varid = -1;
    GET_1D_INT_VAR( "node_num_map", varid, ptr );
    if( -1 == varid )
    {
        MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting node number map data" );
    }

    // Read in the deformations.
    std::vector< std::vector< double > > deformed_arrays( 3 );
    std::vector< std::vector< double > > orig_coords( 3 );
    deformed_arrays[0].reserve( numberNodes_loading );
    deformed_arrays[1].reserve( numberNodes_loading );
    deformed_arrays[2].reserve( numberNodes_loading );
    orig_coords[0].reserve( numberNodes_loading );
    orig_coords[1].reserve( numberNodes_loading );
    orig_coords[2].reserve( numberNodes_loading );
    size_t start[2] = { static_cast< size_t >( time_step - 1 ), 0 },
           count[2] = { 1, static_cast< size_t >( numberNodes_loading ) };
    std::vector< int > dims;
    int coordx = -1, coordy = -1, coordz = -1;
    GET_VAR( "vals_nod_var1", coordx, dims );
    GET_VAR( "vals_nod_var2", coordy, dims );
    if( numberDimensions_loading == 3 ) GET_VAR( "vals_nod_var3", coordz, dims );
    if( -1 == coordx || -1 == coordy || ( numberDimensions_loading == 3 && -1 == coordz ) )
    {
        MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting coords variable" );
    }

    fail = nc_get_vara_double( ncFile, coordx, start, count, &deformed_arrays[0][0] );
    if( NC_NOERR != fail )
    {
        MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting x deformation array" );
    }

    fail = nc_get_vara_double( ncFile, coordy, start, count, &deformed_arrays[1][0] );
    if( NC_NOERR != fail )
    {
        MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting y deformation array" );
    }
    if( numberDimensions_loading == 3 )
    {
        fail = nc_get_vara_double( ncFile, coordz, start, count, &deformed_arrays[2][0] );
        if( NC_NOERR != fail )
        {
            MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting z deformation array" );
        }
    }

    int coord1 = -1, coord2 = -1, coord3 = -1;
    GET_1D_DBL_VAR( "coordx", coord1, orig_coords[0] );
    if( -1 == coord1 )
    {
        MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting x coord array" );
    }
    GET_1D_DBL_VAR( "coordy", coord2, orig_coords[1] );
    if( -1 == coord2 )
    {
        MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting y coord array" );
    }
    if( numberDimensions_loading == 3 )
    {
        GET_1D_DBL_VAR( "coordz", coord3, orig_coords[2] );
        if( -1 == coord3 )
        {
            MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting z coord array" );
        }
    }

    // b. Deal with DB file : get node info. according to node_num_map.
    if( tokens[0] != "coord" && tokens[0] != "COORD" ) return MB_NOT_IMPLEMENTED;

    if( strcmp( op, "set" ) != 0 && strcmp( op, " set" ) != 0 ) return MB_NOT_IMPLEMENTED;

    // Two methods of matching nodes (id vs. proximity)
    const bool match_node_ids = true;

    // Get nodes in cubit file
    Range cub_verts;
    rval = mdbImpl->get_entities_by_type( cub_file_set, MBVERTEX, cub_verts );MB_CHK_ERR( rval );
    std::cout << "  cub_file_set contains " << cub_verts.size() << " nodes." << std::endl;

    // Some accounting
    std::cout << "  exodus file contains " << numberNodes_loading << " nodes." << std::endl;
    double max_magnitude     = 0;
    double average_magnitude = 0;
    int found                = 0;
    int lost                 = 0;
    std::map< int, EntityHandle > cub_verts_id_map;
    AdaptiveKDTree kdtree( mdbImpl );
    EntityHandle root;

    // Should not use cub verts unless they have been matched. Place in a map
    // for fast handle_by_id lookup.
    std::map< int, EntityHandle > matched_cub_vert_id_map;

    // Place cub verts in a map for searching by id
    if( match_node_ids )
    {
        std::vector< int > cub_ids( cub_verts.size() );
        rval = mdbImpl->tag_get_data( mGlobalIdTag, cub_verts, &cub_ids[0] );MB_CHK_ERR( rval );
        for( unsigned i = 0; i != cub_verts.size(); ++i )
        {
            cub_verts_id_map.insert( std::pair< int, EntityHandle >( cub_ids[i], cub_verts[i] ) );
        }

        // Place cub verts in a kdtree for searching by proximity
    }
    else
    {
        FileOptions tree_opts( "MAX_PER_LEAF=1;SPLITS_PER_DIR=1;CANDIDATE_PLANE_SET=0" );
        rval = kdtree.build_tree( cub_verts, &root, &tree_opts );MB_CHK_ERR( rval );
        AdaptiveKDTreeIter tree_iter;
        rval = kdtree.get_tree_iterator( root, tree_iter );MB_CHK_ERR( rval );
    }

    // For each exo vert, find the matching cub vert
    for( int i = 0; i < numberNodes_loading; ++i )
    {
        int exo_id = ptr[i];
        CartVect exo_coords( orig_coords[0][i], orig_coords[1][i], orig_coords[2][i] );
        EntityHandle cub_vert = 0;
        bool found_match      = false;

        // By id
        if( match_node_ids )
        {
            std::map< int, EntityHandle >::iterator i_iter;
            i_iter = cub_verts_id_map.find( exo_id );
            if( i_iter != cub_verts_id_map.end() )
            {
                found_match = true;
                cub_vert    = i_iter->second;
            }

            // By proximity
        }
        else
        {
            // The MAX_NODE_DIST is the farthest distance to  search for a node.
            // For the 1/12th symmetry 85 pin model, the max node dist could not be less
            // than 1e-1 (March 26, 2010).
            const double MAX_NODE_DIST = 1e-1;

            std::vector< EntityHandle > leaves;
            double min_dist = MAX_NODE_DIST;
            rval            = kdtree.distance_search( exo_coords.array(), MAX_NODE_DIST, leaves );MB_CHK_ERR( rval );
            for( std::vector< EntityHandle >::const_iterator j = leaves.begin(); j != leaves.end(); ++j )
            {
                std::vector< EntityHandle > leaf_verts;
                rval = mdbImpl->get_entities_by_type( *j, MBVERTEX, leaf_verts );MB_CHK_ERR( rval );
                for( std::vector< EntityHandle >::const_iterator k = leaf_verts.begin(); k != leaf_verts.end(); ++k )
                {
                    CartVect orig_cub_coords, difference;
                    rval = mdbImpl->get_coords( &( *k ), 1, orig_cub_coords.array() );MB_CHK_ERR( rval );
                    difference  = orig_cub_coords - exo_coords;
                    double dist = difference.length();
                    if( dist < min_dist )
                    {
                        min_dist = dist;
                        cub_vert = *k;
                    }
                }
            }
            if( 0 != cub_vert ) found_match = true;
        }

        // If a match is found, update it with the deformed coords from the exo file.
        if( found_match )
        {
            CartVect updated_exo_coords;
            matched_cub_vert_id_map.insert( std::pair< int, EntityHandle >( exo_id, cub_vert ) );
            updated_exo_coords[0] = orig_coords[0][i] + deformed_arrays[0][i];
            updated_exo_coords[1] = orig_coords[1][i] + deformed_arrays[1][i];

            if( numberDimensions_loading == 3 ) updated_exo_coords[2] = orig_coords[2][i] + deformed_arrays[2][i];
            rval = mdbImpl->set_coords( &cub_vert, 1, updated_exo_coords.array() );MB_CHK_ERR( rval );
            ++found;

            double magnitude =
                sqrt( deformed_arrays[0][i] * deformed_arrays[0][i] + deformed_arrays[1][i] * deformed_arrays[1][i] +
                      deformed_arrays[2][i] * deformed_arrays[2][i] );
            if( magnitude > max_magnitude ) max_magnitude = magnitude;
            average_magnitude += magnitude;
        }
        else
        {
            ++lost;
            std::cout << "cannot match exo node " << exo_id << " " << exo_coords << std::endl;
        }
    }

    // Exo verts that could not be matched have already been printed. Now print the
    // cub verts that could not be matched.
    if( matched_cub_vert_id_map.size() < cub_verts.size() )
    {
        Range unmatched_cub_verts = cub_verts;
        for( std::map< int, EntityHandle >::const_iterator i = matched_cub_vert_id_map.begin();
             i != matched_cub_vert_id_map.end(); ++i )
        {
            unmatched_cub_verts.erase( i->second );
        }

        for( Range::const_iterator i = unmatched_cub_verts.begin(); i != unmatched_cub_verts.end(); ++i )
        {
            int cub_id;
            rval = mdbImpl->tag_get_data( mGlobalIdTag, &( *i ), 1, &cub_id );MB_CHK_ERR( rval );

            CartVect cub_coords;
            rval = mdbImpl->get_coords( &( *i ), 1, cub_coords.array() );MB_CHK_ERR( rval );
            std::cout << "cannot match cub node " << cub_id << " " << cub_coords << std::endl;
        }

        std::cout << "  " << unmatched_cub_verts.size() << " nodes from the cub file could not be matched."
                  << std::endl;
    }

    // Summarize statistics
    std::cout << "  " << found << " nodes from the exodus file were matched in the cub_file_set ";
    if( match_node_ids )
    {
        std::cout << "by id." << std::endl;
    }
    else
    {
        std::cout << "by proximity." << std::endl;
    }

    // Fail if all of the nodes could not be matched.
    if( 0 != lost )
    {
        std::cout << "Error:  " << lost << " nodes from the exodus file could not be matched." << std::endl;
        // return MB_FAILURE;
    }
    std::cout << "  maximum node displacement magnitude: " << max_magnitude << " cm" << std::endl;
    std::cout << "  average node displacement magnitude: " << average_magnitude / found << " cm" << std::endl;

    // *******************************************************************
    // Remove dead elements from the MOAB instance.
    // *******************************************************************

    // How many element variables are in the file?
    int n_elem_var;
    GET_DIM( ncdim, "num_elem_var", n_elem_var );

    // Get element variable names
    varid       = -1;
    int cstatus = nc_inq_varid( ncFile, "name_elem_var", &varid );
    std::vector< char > names_memory( n_elem_var * max_str_length );
    std::vector< char* > names( n_elem_var );
    for( int i = 0; i < n_elem_var; ++i )
        names[i] = &names_memory[i * max_str_length];
    if( cstatus != NC_NOERR || varid == -1 )
    {
        std::cout << "ReadNCDF: name_elem_var does not exist" << std::endl;
        return MB_FAILURE;
    }
    int status = nc_get_var_text( ncFile, varid, &names_memory[0] );
    if( NC_NOERR != status )
    {
        MB_SET_ERR( MB_FAILURE, "ReadNCDF: Problem getting element variable names" );
    }

    // Is one of the element variable names "death_status"? If so, get its index
    // in the element variable array.
    int death_index;
    bool found_death_index = false;
    for( int i = 0; i < n_elem_var; ++i )
    {
        std::string temp( names[i] );
        std::string::size_type pos0 = temp.find( "death_status" );
        std::string::size_type pos1 = temp.find( "Death_Status" );
        std::string::size_type pos2 = temp.find( "DEATH_STATUS" );
        if( std::string::npos != pos0 || std::string::npos != pos1 || std::string::npos != pos2 )
        {
            found_death_index = true;
            death_index       = i + 1;  // NetCDF variables start with 1
            break;
        }
    }
    if( !found_death_index )
    {
        MB_SET_ERR( MB_FAILURE, "ReadNCDF: Problem getting index of death_status variable" );
    }

    // The exodus header has already been read. This contains the number of element
    // blocks.

    // Dead elements are listed by block. Read the block headers to determine how
    // many elements are in each block.
    rval = read_block_headers( blocks_to_load, num_blocks );
    if( MB_FAILURE == rval )
    {
        MB_SET_ERR( MB_FAILURE, "ReadNCDF: Problem reading block headers" );
    }

    // Dead elements from the Exodus file can be located in the cub_file_set by id
    // or by connectivity. Currently, finding elements by id requires careful book
    // keeping when constructing the model in Cubit. To avoid this, one can match
    // elements by connectivity instead.
    const bool match_elems_by_connectivity = true;

    // Get the element id map. The ids in the map are from the elements in the blocks.
    // elem_num_map(blk1 elem ids, blk2 elem ids, blk3 elem ids, ...)
    std::vector< int > elem_ids( numberNodes_loading );
    if( !match_elems_by_connectivity )
    {
        GET_1D_INT_VAR( "elem_num_map", varid, elem_ids );
        if( -1 == varid )
        {
            MB_SET_ERR( MB_FAILURE, "ReadNCDF: Problem getting element number map data" );
        }
    }

    // For each block
    int first_elem_id_in_block = 0;
    int block_count            = 1;  // NetCDF variables start with 1
    int total_elems            = 0;
    int total_dead_elems       = 0;
    for( std::vector< ReadBlockData >::iterator i = blocksLoading.begin(); i != blocksLoading.end(); ++i )
    {
        // Get the ncdf connect variable
        std::string temp_string( "connect" );
        std::stringstream temp_ss;
        temp_ss << block_count;
        temp_string += temp_ss.str();
        temp_string += "\0";
        std::vector< int > ldims;
        GET_VAR( temp_string.c_str(), nc_var, ldims );
        if( !nc_var )
        {
            MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting connect variable" );
        }
        // The element type is an attribute of the connectivity variable
        nc_type att_type;
        size_t att_len;
        fail = nc_inq_att( ncFile, nc_var, "elem_type", &att_type, &att_len );
        if( NC_NOERR != fail )
        {
            MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting elem type attribute" );
        }
        std::vector< char > dum_str( att_len + 1 );
        fail = nc_get_att_text( ncFile, nc_var, "elem_type", &dum_str[0] );
        if( NC_NOERR != fail )
        {
            MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting elem type" );
        }
        ExoIIElementType elem_type = ExoIIUtil::static_element_name_to_type( &dum_str[0] );
        ( *i ).elemType            = elem_type;
        const EntityType mb_type   = ExoIIUtil::ExoIIElementMBEntity[( *i ).elemType];

        // Get the number of nodes per element
        unsigned int nodes_per_element = ExoIIUtil::VerticesPerElement[( *i ).elemType];

        // Read the connectivity into that memory.
        // int exo_conn[i->numElements][nodes_per_element];
        int* exo_conn = new int[i->numElements * nodes_per_element];
        // NcBool status = temp_var->get(&exo_conn[0][0], i->numElements, nodes_per_element);
        size_t lstart[2] = { 0, 0 }, lcount[2];
        lcount[0]        = i->numElements;
        lcount[1]        = nodes_per_element;
        fail             = nc_get_vara_int( ncFile, nc_var, lstart, lcount, exo_conn );
        if( NC_NOERR != fail )
        {
            delete[] exo_conn;
            MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting connectivity" );
        }

        // Get the death_status at the correct time step.
        std::vector< double > death_status( i->numElements );  // It seems wrong, but it uses doubles
        std::string array_name( "vals_elem_var" );
        temp_ss.str( "" );  // stringstream won't clear by temp.clear()
        temp_ss << death_index;
        array_name += temp_ss.str();
        array_name += "eb";
        temp_ss.str( "" );  // stringstream won't clear by temp.clear()
        temp_ss << block_count;
        array_name += temp_ss.str();
        array_name += "\0";
        GET_VAR( array_name.c_str(), nc_var, ldims );
        if( !nc_var )
        {
            MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting death status variable" );
        }
        lstart[0] = time_step - 1;
        lstart[1] = 0;
        lcount[0] = 1;
        lcount[1] = i->numElements;
        status    = nc_get_vara_double( ncFile, nc_var, lstart, lcount, &death_status[0] );
        if( NC_NOERR != status )
        {
            delete[] exo_conn;
            MB_SET_ERR( MB_FAILURE, "ReadNCDF: Problem setting time step for death_status" );
        }

        // Look for dead elements. If there is too many dead elements and this starts
        // to take too long, I should put the elems in a kd-tree for more efficient
        // searching. Alternatively I could get the exo connectivity and match nodes.
        int dead_elem_counter = 0, missing_elem_counter = 0;
        for( int j = 0; j < i->numElements; ++j )
        {
            if( 1 != death_status[j] )
            {
                Range cub_elem, cub_nodes;
                if( match_elems_by_connectivity )
                {
                    // Get exodus nodes for the element
                    std::vector< int > elem_conn( nodes_per_element );
                    for( unsigned int k = 0; k < nodes_per_element; ++k )
                    {
                        // elem_conn[k] = exo_conn[j][k];
                        elem_conn[k] = exo_conn[j * nodes_per_element + k];
                    }
                    // Get the ids of the nodes (assume we are matching by id)
                    // Remember that the exodus array locations start with 1 (not 0).
                    std::vector< int > elem_conn_node_ids( nodes_per_element );
                    for( unsigned int k = 0; k < nodes_per_element; ++k )
                    {
                        elem_conn_node_ids[k] = ptr[elem_conn[k] - 1];
                    }
                    // Get the cub_file_set nodes by id
                    // The map is a log search and takes almost no time.
                    // MOAB's linear tag search takes 5-10 minutes.
                    for( unsigned int k = 0; k < nodes_per_element; ++k )
                    {
                        std::map< int, EntityHandle >::iterator k_iter;
                        k_iter = matched_cub_vert_id_map.find( elem_conn_node_ids[k] );

                        if( k_iter == matched_cub_vert_id_map.end() )
                        {
                            std::cout << "ReadNCDF: Found no cub node with id=" << elem_conn_node_ids[k]
                                      << ", but expected to find only 1." << std::endl;
                            break;
                        }
                        cub_nodes.insert( k_iter->second );
                    }

                    if( nodes_per_element != cub_nodes.size() )
                    {
                        std::cout << "ReadNCDF: nodes_per_elemenet != cub_nodes.size()" << std::endl;
                        delete[] exo_conn;
                        return MB_INVALID_SIZE;
                    }

                    // Get the cub_file_set element with the same nodes
                    int to_dim = CN::Dimension( mb_type );
                    rval       = mdbImpl->get_adjacencies( cub_nodes, to_dim, false, cub_elem );
                    if( MB_SUCCESS != rval )
                    {
                        std::cout << "ReadNCDF: could not get dead cub element" << std::endl;
                        delete[] exo_conn;
                        return rval;
                    }

                    // Pronto/Presto renumbers elements, so matching cub and exo elements by
                    // id is not possible at this time.
                }
                else
                {
                    // Get dead element's id
                    int elem_id = elem_ids[first_elem_id_in_block + j];
                    void* id[]  = { &elem_id };
                    // Get the element by id
                    rval = mdbImpl->get_entities_by_type_and_tag( cub_file_set, mb_type, &mGlobalIdTag, id, 1, cub_elem,
                                                                  Interface::INTERSECT );
                    if( MB_SUCCESS != rval )
                    {
                        delete[] exo_conn;
                        return rval;
                    }
                }

                if( 1 == cub_elem.size() )
                {
                    // Delete the dead element from the cub file. It will be removed from sets
                    // ONLY if they are tracking meshsets.
                    rval = mdbImpl->remove_entities( cub_file_set, cub_elem );
                    if( MB_SUCCESS != rval )
                    {
                        delete[] exo_conn;
                        return rval;
                    }
                    rval = mdbImpl->delete_entities( cub_elem );
                    if( MB_SUCCESS != rval )
                    {
                        delete[] exo_conn;
                        return rval;
                    }
                }
                else
                {
                    std::cout << "ReadNCDF: Should have found 1 element with  type=" << mb_type
                              << " in cub_file_set, but instead found " << cub_elem.size() << std::endl;
                    rval = mdbImpl->list_entities( cub_nodes );
                    ++missing_elem_counter;
                    delete[] exo_conn;
                    return MB_FAILURE;
                }
                ++dead_elem_counter;
            }
        }
        // Print some statistics
        temp_ss.str( "" );
        temp_ss << i->blockId;
        total_dead_elems += dead_elem_counter;
        total_elems += i->numElements;
        std::cout << "  Block " << temp_ss.str() << " has " << dead_elem_counter << "/" << i->numElements
                  << " dead elements." << std::endl;
        if( 0 != missing_elem_counter )
        {
            std::cout << "    " << missing_elem_counter
                      << " dead elements in this block were not found in the cub_file_set." << std::endl;
        }

        // Advance the pointers into element ids and block_count. Memory cleanup.
        first_elem_id_in_block += i->numElements;
        ++block_count;
        delete[] exo_conn;
    }

    std::cout << " Total: " << total_dead_elems << "/" << total_elems << " dead elements." << std::endl;

    // Close the file
    fail = nc_close( ncFile );
    if( NC_NOERR != fail )
    {
        MB_SET_ERR( MB_FAILURE, "Trouble closing file" );
    }

    ncFile = 0;
    return MB_SUCCESS;
}

void ReadNCDF::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 );
    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();
    }
}

}  // namespace moab