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