Output Exodus File for VERDE. More...

#include <ReadNCDF.hpp>

Inheritance diagram for moab::ReadNCDF:

Public Member Functions
void	tokenize (const std::string &str, std::vector< std::string > &tokens, const char *delimiters)
ErrorCode	load_file (const char file_name, const EntityHandle file_set, const FileOptions &opts, const SubsetList subset_list=0, const Tag file_id_tag=0)
	load an ExoII file
ErrorCode	read_tag_values (const char file_name, const char tag_name, const FileOptions &opts, std::vector< int > &tag_values_out, const SubsetList *subset_list=0)
	Read tag values from a file.
	ReadNCDF (Interface *impl=NULL)
	Constructor.
virtual	~ReadNCDF ()
	Destructor.
ErrorCode	update (const char exodus_file_name, const FileOptions &opts, const int num_blocks, const int blocks_to_load, const EntityHandle file_set)
Static Public Member Functions
static ReaderIface *	factory (Interface *)
Private Member Functions
bool	dimension_exists (const char *attrib_name)
void	reset ()
ErrorCode	read_exodus_header ()
	read the header from the ExoII file
ErrorCode	read_nodes (const Tag *file_id_tag)
	read the nodes
ErrorCode	read_block_headers (const int *blocks_to_load, const int num_blocks)
	read block headers, containing info about element type, number, etc.
ErrorCode	read_elements (const Tag *file_id_tag)
	read the element blocks
ErrorCode	read_global_ids ()
	read in the global element ids
ErrorCode	read_nodesets ()
	read the nodesets into meshsets
ErrorCode	read_sidesets ()
	read the sidesets (does nothing for now)
int	exodus_file ()
	exodus file bound to this object
int	number_dimensions ()
	number of dimensions in this exo file
ErrorCode	get_type (char *exo_element_type, EntityType &elem_type)
	map a character exodusII element type to a TSTT type & topology
ErrorCode	get_type (EntityType &elem_type, std::string &exo_element_type)
ErrorCode	read_qa_records (EntityHandle file_set)
ErrorCode	read_qa_information (std::vector< std::string > &qa_record_list)
ErrorCode	read_qa_string (char *string, int record_number, int record_position)
ErrorCode	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)
ErrorCode	find_side_element_type (const int element_id, ExoIIElementType &type, ReadBlockData &block_data, int &df_index, int side_id)
ErrorCode	create_sideset_element (const std::vector< EntityHandle > &, EntityType, EntityHandle &, int &)
	creates an element with the given connectivity
int	get_number_nodes (EntityHandle handle)
Private Attributes
ReadUtilIface *	readMeshIface
Interface *	mdbImpl
	interface instance
int	ncFile
int	CPU_WORD_SIZE
int	IO_WORD_SIZE
int	vertexOffset
	int to offset vertex ids with
std::string	exodusFile
	file name
int	numberNodes_loading
	number of nodes in the current exoII file
int	numberElements_loading
	number of elements in the current exoII file
int	numberElementBlocks_loading
	number of blocks in the current exoII file
int	numberNodeSets_loading
	number of nodesets in the current exoII file
int	numberSideSets_loading
	number of sidesets in the current exoII file
int	numberDimensions_loading
EntityHandle	mCurrentMeshHandle
	Meshset Handle for the mesh that is currently being read.
std::vector< char >	nodesInLoadedBlocks
std::vector< ReadBlockData >	blocksLoading
Tag	mMaterialSetTag
Tag	mDirichletSetTag
Tag	mNeumannSetTag
Tag	mHasMidNodesTag
Tag	mDistFactorTag
Tag	mGlobalIdTag
Tag	mQaRecordTag
int	max_line_length
int	max_str_length
Range	initRange
	range of entities in initial mesh, before this read

Detailed Description

Output Exodus File for VERDE.

Definition at line 60 of file ReadNCDF.hpp.

Constructor & Destructor Documentation

moab::ReadNCDF::ReadNCDF ( Interface * impl = NULL )

Constructor.

get and cache predefined tag handles

Definition at line 107 of file ReadNCDF.cpp.

    : 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;

  ErrorCode result;
  const int zero = 0, 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);
  result = impl->tag_get_handle(GLOBAL_ID_TAG_NAME, 1, MB_TYPE_INTEGER,
                                mGlobalIdTag, MB_TAG_SPARSE|MB_TAG_CREAT, &zero);
  assert(MB_SUCCESS == result);
#ifdef NDEBUG
  if (MB_SUCCESS == result) {}; // line to avoid compiler warning about unused variable
#endif
  ncFile = 0;
}

moab::ReadNCDF::~ReadNCDF ( ) [virtual]

Destructor.

Definition at line 174 of file ReadNCDF.cpp.

{
  mdbImpl->release_interface(readMeshIface);
}

Member Function Documentation

ErrorCode moab::ReadNCDF::create_sideset_element	(	const std::vector< EntityHandle > &	connectivity,
		EntityType	type,
		EntityHandle &	handle,
		int &	sense
	)		`[private]`

creates an element with the given connectivity

Definition at line 1325 of file ReadNCDF.cpp.

{
  // get adjacent entities
  ErrorCode error = MB_SUCCESS;
  int to_dim = CN::Dimension(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;
    iter = std::find(match_conn.begin(), match_conn.end(), connectivity[0]);
    if(iter == match_conn.end())
      continue;
   
    // rotate to match connectivity 
    std::rotate(match_conn.begin(), iter, match_conn.end());

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

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

      unsigned int k = connectivity.size() - 1;
      for(j=1; j<connectivity.size(); )
      {
        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 == true)
      handle = adj_ent[i];
  }

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

  return error;
}

ErrorCode moab::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
	)		`[private]`

Definition at line 1084 of file ReadNCDF.cpp.

{
  //determine entity type from element id
  int i, k;

  // 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) {
      readMeshIface->report_error("ReadNCDF:: Problem getting dist fact variable.");
      return MB_FAILURE;
    }
  }

  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;
  for(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 (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(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 (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(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(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(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 (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(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 (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(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(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 (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(k=0; k<2; k++)
            dist_factor_vector.push_back( temp_dist_factor_vector[df_index++] );
        }
      }
    }
  }

  return MB_SUCCESS; 
}

bool moab::ReadNCDF::dimension_exists ( const char * attrib_name ) [private]

int moab::ReadNCDF::exodus_file ( ) [private]

exodus file bound to this object

ReaderIface * moab::ReadNCDF::factory ( Interface * iface ) [static]

Definition at line 104 of file ReadNCDF.cpp.

  { return new ReadNCDF( iface ); }

ErrorCode moab::ReadNCDF::find_side_element_type	(	const int	element_id,
		ExoIIElementType &	type,
		ReadBlockData &	block_data,
		int &	df_index,
		int	side_id
	)		`[private]`

Definition at line 1403 of file ReadNCDF.cpp.

{
  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;
}

int moab::ReadNCDF::get_number_nodes ( EntityHandle handle ) [private]

ErrorCode moab::ReadNCDF::get_type	(	char *	exo_element_type,
		EntityType &	elem_type
	)		`[private]`

map a character exodusII element type to a TSTT type & topology

ErrorCode moab::ReadNCDF::get_type	(	EntityType &	elem_type,
		std::string &	exo_element_type
	)		`[private]`

ErrorCode moab::ReadNCDF::load_file	(	const char *	file_name,
		const EntityHandle *	file_set,
		const FileOptions &	opts,
		const SubsetList *	subset_list = `0`,
		const Tag *	file_id_tag = `0`
	)		`[virtual]`

load an ExoII file

Implements moab::ReaderIface.

Definition at line 237 of file ReadNCDF.cpp.

{
  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) ) {
      readMeshIface->report_error( "ExodusII reader supports subset read only by material ID." );
      return MB_UNSUPPORTED_OPERATION;
    }
    if (subset_list->num_parts) {
      readMeshIface->report_error( "ExodusII reader does not support mesh partitioning");
      return MB_UNSUPPORTED_OPERATION;
    }
    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)
  {
    readMeshIface->report_error("ReadNCDF:: problem opening Netcdf/Exodus II file %s",exodus_file_name);
    return MB_FILE_DOES_NOT_EXIST;
  }

    // 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. 
  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)
  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???

  ncFile = 0;
  return MB_SUCCESS;
}

int moab::ReadNCDF::number_dimensions ( ) [inline, private]

number of dimensions in this exo file

Definition at line 232 of file ReadNCDF.hpp.

{
   return numberDimensions_loading;
}

ErrorCode moab::ReadNCDF::read_block_headers	(	const int *	blocks_to_load,
		const int	num_blocks
	)		`[private]`

read block headers, containing info about element type, number, etc.

Definition at line 483 of file ReadNCDF.cpp.

{
     
    // 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) {
    readMeshIface->report_error("ReadNCDF:: Problem getting eb_prop1 variable.");
    return MB_FAILURE;
  }

  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) 
        != block_ids.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 moab::ReadNCDF::read_elements ( const Tag * file_id_tag ) [private]

read the element blocks

Definition at line 558 of file ReadNCDF.cpp.

{
    // read in elements
  
  int result = 0;

  // intialize the nodeInLoadedBlocks vector
  nodesInLoadedBlocks.resize(numberNodes_loading+1);
  memset(&nodesInLoadedBlocks[0], 0, (numberNodes_loading+1)*sizeof(char));

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


  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 (!nc_var) {
      readMeshIface->report_error("ReadNCDF:: Problem getting connect variable.");
      return MB_FAILURE;
    }
    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) {
      readMeshIface->report_error("ReadNCDF:: Problem getting elem type attribute.");
      return MB_FAILURE;
    }
    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) {
      readMeshIface->report_error("ReadNCDF:: Problem getting elem type.");
      return MB_FAILURE;
    }
    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];

    // 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 connetivity 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) {
      readMeshIface->report_error("ReadNCDF:: Problem getting connectivity.");
      return MB_FAILURE;
    }

      // 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()) {
        readMeshIface->report_error( "Invalid node ID in block connectivity\n");
        return MB_FAILURE;
      }
      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 )
    {
      readMeshIface->report_error("ReadNCDF:: error getting element connectivity for block %i",
          block_id);
      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;

    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 );
    }
  }

  return MB_SUCCESS;
}

ErrorCode moab::ReadNCDF::read_exodus_header ( ) [private]

read the header from the ExoII file

Definition at line 327 of file ReadNCDF.cpp.

{
  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) {
    readMeshIface->report_error("ReadNCDF: File contains %d dims but NetCDF library supports only %d\n",
                                ndims, (int)NC_MAX_DIMS);
    return MB_FAILURE;
  }
  int nvars;
  fail = nc_inq_nvars(ncFile, &nvars);
  if (nvars > NC_MAX_VARS) {
    readMeshIface->report_error("ReadNCDF: File contains %d vars but NetCDF library supports only %d\n",
                                nvars, (int)NC_MAX_VARS);
    return MB_FAILURE;
  }
  
    // 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) {
    readMeshIface->report_error("ReadNCDF:: Problem getting floating_point_word_size attribute.");
    return MB_FAILURE;
  }
  if (att_type != NC_INT || att_len != 1) {
    readMeshIface->report_error("ReadNCDF:: Word size didn't have type int or size 1.");
    return MB_FAILURE;
  }
  fail = nc_get_att_int(ncFile, NC_GLOBAL, "floating_point_word_size", &IO_WORD_SIZE);
  if (NC_NOERR != fail) {
    readMeshIface->report_error("ReadNCDF:: Trouble getting word size.");
    return MB_FAILURE;
  }

    // exodus version
  fail = nc_inq_att(ncFile, NC_GLOBAL, "version", &att_type, &att_len);
  if (NC_NOERR != fail) {
    readMeshIface->report_error("ReadNCDF:: Problem getting version attribute.");
    return MB_FAILURE;
  }
  if (att_type != NC_FLOAT || att_len != 1) {
    readMeshIface->report_error("ReadNCDF:: Version didn't have type float or size 1.");
    return MB_FAILURE;
  }
    // 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_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) {
    readMeshIface->report_error("ReadNCDF:: Problem getting title attribute.");
    return MB_FAILURE;
  }
  if (att_type != NC_CHAR) {
    readMeshIface->report_error("ReadNCDF:: title didn't have type char.");
    return MB_FAILURE;
  }
  char *title = new char[att_len + 1];
  fail = nc_get_att_text(ncFile, NC_GLOBAL, "title", title);
  if (NC_NOERR != fail) {
    readMeshIface->report_error("ReadNCDF:: trouble getting title.");
    delete[] title;
    return MB_FAILURE;
  }
  delete[] title;

  return MB_SUCCESS;
}

ErrorCode moab::ReadNCDF::read_global_ids ( ) [private]

read in the global element ids

Definition at line 708 of file ReadNCDF.cpp.

{
    // 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->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;
        }
      }
    }
  }

    // 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)
      readMeshIface->report_error("ReadNCDF:: Problem setting node global ids.");
  }
  
  return MB_SUCCESS;
}

ErrorCode moab::ReadNCDF::read_nodes ( const Tag * file_id_tag ) [private]

read the nodes

Definition at line 416 of file ReadNCDF.cpp.

{
    // 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) {
        readMeshIface->report_error("ReadNCDF:: Problem getting %c coord array.", 'x'+d);
        return MB_FAILURE;
      }
    }
  }
  // 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) {
        readMeshIface->report_error("ReadNCDF:: Problem getting %c coord variable.", 'x'+d);
        return MB_FAILURE;
      }
      
      fail = nc_get_vara_double(ncFile, coord, start, &count[1], arrays[d]);
      if (NC_NOERR != fail) {
        readMeshIface->report_error("ReadNCDF:: Problem getting %c coord array.", 'x'+d);
        return MB_FAILURE;
      }
    }
  }
  
  // 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 moab::ReadNCDF::read_nodesets ( ) [private]

read the nodesets into meshsets

Definition at line 754 of file ReadNCDF.cpp.

{
    // 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) {
    readMeshIface->report_error("ReadNCDF:: Problem getting ns_prop1 variable.");
    return MB_FAILURE;
  }

    // 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) {
        readMeshIface->report_error("ReadNCDF:: Problem getting dist fact variable.");
        return MB_FAILURE;
      }
    }

      // 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) {
      readMeshIface->report_error("ReadNCDF:: Problem getting nodeset node variable.");
      return MB_FAILURE;
    }

      // 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 moab::ReadNCDF::read_qa_information ( std::vector< std::string > & qa_record_list ) [private]

Definition at line 1474 of file ReadNCDF.cpp.

{
  // 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 moab::ReadNCDF::read_qa_records ( EntityHandle file_set ) [private]

Definition at line 1450 of file ReadNCDF.cpp.

{
  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 moab::ReadNCDF::read_qa_string	(	char *	string,
		int	record_number,
		int	record_position
	)		`[private]`

Definition at line 1499 of file ReadNCDF.cpp.

{
  std::vector<int> dims;
  int temp_var = -1;
  GET_VAR("qa_records", temp_var, dims);
  if (-1 == temp_var) {
    readMeshIface->report_error("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) {
    readMeshIface->report_error("ReadNCDF:: Problem setting current record number variable position.");
    return MB_FAILURE;
  }
    // get the variable id in the exodus file

  return MB_SUCCESS;
}

ErrorCode moab::ReadNCDF::read_sidesets ( ) [private]

read the sidesets (does nothing for now)

Definition at line 921 of file ReadNCDF.cpp.

{
  // 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) {
    readMeshIface->report_error("ReadNCDF:: Problem getting ss_prop1 variable.");
    return MB_FAILURE;
  }

  // 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) {
      readMeshIface->report_error("ReadNCDF:: Problem getting sideset side variable.");
      return MB_FAILURE;
    }

    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) {
      readMeshIface->report_error("ReadNCDF:: Problem getting sideset elem variable.");
      return MB_FAILURE;
    }

    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[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 moab::ReadNCDF::read_tag_values	(	const char *	file_name,
		const char *	tag_name,
		const FileOptions &	opts,
		std::vector< int > &	tag_values_out,
		const SubsetList *	subset_list = `0`
	)		`[virtual]`

Read tag values from a file.

Read the list if all integer tag values from the file for a tag that is a single integer value per entity.

Parameters:

file_name	The file to read.
tag_name	The tag for which to read values
tag_values_out	Output: The list of tag values.
subset_list	An array of tag name and value sets specifying the subset of the file to read. If multiple tags are specified, the sets that match all tags (intersection) should be read.
subset_list_length	The length of the 'subset_list' array.

Implements moab::ReaderIface.

Definition at line 179 of file ReadNCDF.cpp.

{
  if (subset_list) {
    readMeshIface->report_error( "ExodusII reader supports subset read only by material ID." );
    return MB_UNSUPPORTED_OPERATION;
  }

      // open netcdf/exodus file
  int fail = nc_open(file_name, 0, &ncFile);
  if (NC_NOWRITE != fail) {
    readMeshIface->report_error("ReadNCDF:: problem opening Netcdf/Exodus II file %s",file_name);
    return MB_FILE_DOES_NOT_EXIST;
  }

    // 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) {
      readMeshIface->report_error("Problem getting prop variable.");
      return MB_FAILURE;
    }
  }
  
  return MB_SUCCESS;
}

void moab::ReadNCDF::reset ( ) [private]

Definition at line 154 of file ReadNCDF.cpp.

{
  numberDimensions_loading = -1;
  mCurrentMeshHandle = 0;
  vertexOffset = 0; 

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

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

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

void moab::ReadNCDF::tokenize	(	const std::string &	str,
		std::vector< std::string > &	tokens,
		const char *	delimiters
	)

Definition at line 2139 of file ReadNCDF.cpp.

{
  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();
  }
}

ErrorCode moab::ReadNCDF::update	(	const char *	exodus_file_name,
		const FileOptions &	opts,
		const int	num_blocks,
		const int *	blocks_to_load,
		const EntityHandle	file_set
	)

Definition at line 1526 of file ReadNCDF.cpp.

{
  //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)
  {
    readMeshIface->report_error("ReadNCDF:: Problem reading file options.");
    return MB_FAILURE;
  }
  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")) {
    readMeshIface->report_error("ReadNCDF: invalid operation specified for update");
    return MB_TYPE_OUT_OF_RANGE;
  }
  op = tokens[2].c_str();

      // open netcdf/exodus file
  ncFile = 0;
  int fail = nc_open(exodus_file_name, 0, &ncFile);
  if (!ncFile)
  {
    readMeshIface->report_error("ReadNCDF:: problem opening Netcdf/Exodus II file %s",exodus_file_name);
    return MB_FILE_DOES_NOT_EXIST;
  }

  rval = read_exodus_header();
  if (MB_SUCCESS != rval)
    return 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) {
    readMeshIface->report_error("ReadNCDF:: Problem getting node number map data.");
    return MB_FAILURE;
  }

  // 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)) {
    readMeshIface->report_error("ReadNCDF:: Problem getting coords variable.");
    return MB_FAILURE;
  }

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

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

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

  //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") && strcmp (op, " set"))
    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 );
  if(MB_SUCCESS != rval) return 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] );
    if(MB_SUCCESS != rval) return 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);
    if(MB_SUCCESS != rval) return rval;
    AdaptiveKDTreeIter tree_iter;                                                     
    rval = kdtree.get_tree_iterator( root, tree_iter );
    if(MB_SUCCESS != rval) return 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);
      if(MB_SUCCESS != rval) return 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 );
    if(MB_SUCCESS != rval) return 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() );
          if(MB_SUCCESS != rval) return 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() );
      if(MB_SUCCESS != rval) return 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 );
      if(MB_SUCCESS != rval) return rval;
      CartVect cub_coords;
      rval = mdbImpl->get_coords( &(*i), 1, cub_coords.array() );      
      if(MB_SUCCESS != rval) return 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) {
    readMeshIface->report_error("ReadNCDF: Problem getting element variable names.");
    return MB_FAILURE;
  }

  // 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) {
    std::cout << "ReadNCDF: Problem getting index of death_status variable." << std::endl;
    return MB_FAILURE;
  }
    
  // 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) {
    std::cout << "ReadNCDF: Problem reading block headers." << std::endl;
    return rval;
  }

  // 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 elementsin 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) {
      std::cout << "ReadNCDF: Problem getting element number map data." << std::endl;
      return MB_FAILURE;
    }
  }

  // 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) {
      readMeshIface->report_error("ReadNCDF:: Problem getting connect variable.");
      return MB_FAILURE;
    }
    // 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) {
      readMeshIface->report_error("ReadNCDF:: Problem getting elem type attribute.");
      return MB_FAILURE;
    }
    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) {
      readMeshIface->report_error("ReadNCDF:: Problem getting elem type.");
      return MB_FAILURE;
    }
    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) {
      readMeshIface->report_error("ReadNCDF:: Problem getting connectivity.");
      delete [] exo_conn;
      return MB_FAILURE;
    }

    // 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) {
      readMeshIface->report_error("ReadNCDF:: Problem getting death status variable.");
      return MB_FAILURE;
    }
    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) {
      std::cout << "ReadNCDF: Problem setting time step for death_status." << std::endl;
      delete[] exo_conn;
      return MB_FAILURE;
    }

    // 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;

  return MB_SUCCESS;
}