Mesh Oriented datABase  (version 5.5.0)
An array-based unstructured mesh library
moab::ReadNCDF Class Reference

Output Exodus File for VERDE. More...

#include <ReadNCDF.hpp>

+ Inheritance diagram for moab::ReadNCDF:
+ Collaboration 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 More...
 
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. More...
 
 ReadNCDF (Interface *impl=NULL)
 Constructor. More...
 
virtual ~ReadNCDF ()
 Destructor. More...
 
ErrorCode update (const char *exodus_file_name, const FileOptions &opts, const int num_blocks, const int *blocks_to_load, const EntityHandle file_set)
 
- Public Member Functions inherited from moab::ReaderIface
virtual ~ReaderIface ()
 

Static Public Member Functions

static ReaderIfacefactory (Interface *)
 

Private Member Functions

bool dimension_exists (const char *attrib_name)
 
void reset ()
 
ErrorCode read_exodus_header ()
 read the header from the ExoII file More...
 
ErrorCode read_nodes (const Tag *file_id_tag)
 read the nodes More...
 
ErrorCode read_face_blocks_headers ()
 
ErrorCode read_block_headers (const int *blocks_to_load, const int num_blocks)
 read block headers, containing info about element type, number, etc. More...
 
ErrorCode read_polyhedra_faces ()
 
ErrorCode read_elements (const Tag *file_id_tag)
 read the element blocks More...
 
ErrorCode read_global_ids ()
 read in the global element ids More...
 
ErrorCode read_nodesets ()
 read the nodesets into meshsets More...
 
ErrorCode read_sidesets ()
 read the sidesets (does nothing for now) More...
 
int exodus_file ()
 exodus file bound to this object More...
 
int number_dimensions ()
 number of dimensions in this exo file More...
 
ErrorCode get_type (char *exo_element_type, EntityType &elem_type)
 map a character exodusII element type to a TSTT type & topology More...
 
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 More...
 
int get_number_nodes (EntityHandle handle)
 

Private Attributes

ReadUtilIfacereadMeshIface
 
InterfacemdbImpl
 interface instance More...
 
int ncFile
 
int CPU_WORD_SIZE
 
int IO_WORD_SIZE
 
int vertexOffset
 int to offset vertex ids with More...
 
std::string exodusFile
 file name More...
 
int numberNodes_loading
 number of nodes in the current exoII file More...
 
int numberElements_loading
 number of elements in the current exoII file More...
 
int numberElementBlocks_loading
 number of blocks in the current exoII file More...
 
int numberFaceBlocks_loading
 number of face blocks in the current exoII file (used for polyhedra) More...
 
int numberNodeSets_loading
 number of nodesets in the current exoII file More...
 
int numberSideSets_loading
 number of sidesets in the current exoII file More...
 
int numberDimensions_loading
 
EntityHandle mCurrentMeshHandle
 Meshset Handle for the mesh that is currently being read. More...
 
std::vector< char > nodesInLoadedBlocks
 
std::vector< ReadBlockDatablocksLoading
 
std::vector< EntityHandlepolyfaces
 
std::vector< ReadFaceBlockDatafaceBlocksLoading
 
Tag mMaterialSetTag
 Cached tags for reading. Note that all these tags are defined when the core is initialized. More...
 
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 More...
 

Detailed Description

Output Exodus File for VERDE.

Definition at line 73 of file ReadNCDF.hpp.

Constructor & Destructor Documentation

◆ ReadNCDF()

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

Constructor.

Get and cache predefined tag handles

Definition at line 140 of file ReadNCDF.cpp.

140  : mdbImpl( impl ), max_line_length( -1 ), max_str_length( -1 )
141 {
142  assert( impl != NULL );
143  reset();
144 
145  impl->query_interface( readMeshIface );
146 
147  // Initialize in case tag_get_handle fails below
148  mMaterialSetTag = 0;
149  mDirichletSetTag = 0;
150  mNeumannSetTag = 0;
151  mHasMidNodesTag = 0;
152  mDistFactorTag = 0;
153  mQaRecordTag = 0;
154  mGlobalIdTag = 0;
155 
156  //! Get and cache predefined tag handles
157  ErrorCode result;
158  const int negone = -1;
159  result = impl->tag_get_handle( MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mMaterialSetTag,
160  MB_TAG_SPARSE | MB_TAG_CREAT, &negone );
161  assert( MB_SUCCESS == result );
162  result = impl->tag_get_handle( DIRICHLET_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mDirichletSetTag,
163  MB_TAG_SPARSE | MB_TAG_CREAT, &negone );
164  assert( MB_SUCCESS == result );
165  result = impl->tag_get_handle( NEUMANN_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mNeumannSetTag,
166  MB_TAG_SPARSE | MB_TAG_CREAT, &negone );
167  assert( MB_SUCCESS == result );
168  const int mids[] = { -1, -1, -1, -1 };
169  result = impl->tag_get_handle( HAS_MID_NODES_TAG_NAME, 4, MB_TYPE_INTEGER, mHasMidNodesTag,
170  MB_TAG_SPARSE | MB_TAG_CREAT, mids );
171  assert( MB_SUCCESS == result );
172  result = impl->tag_get_handle( "distFactor", 0, MB_TYPE_DOUBLE, mDistFactorTag,
174  assert( MB_SUCCESS == result );
175  result = impl->tag_get_handle( "qaRecord", 0, MB_TYPE_OPAQUE, mQaRecordTag,
177  assert( MB_SUCCESS == result );
178  mGlobalIdTag = impl->globalId_tag();
179 
180  assert( MB_SUCCESS == result );
181 #ifdef NDEBUG
182  if( MB_SUCCESS == result )
183  {
184  }; // Line to avoid compiler warning about unused variable
185 #endif
186  ncFile = 0;
187 }

References DIRICHLET_SET_TAG_NAME, ErrorCode, moab::Interface::globalId_tag(), HAS_MID_NODES_TAG_NAME, MATERIAL_SET_TAG_NAME, MB_SUCCESS, MB_TAG_CREAT, MB_TAG_SPARSE, MB_TAG_VARLEN, MB_TYPE_DOUBLE, MB_TYPE_INTEGER, MB_TYPE_OPAQUE, mDirichletSetTag, mDistFactorTag, mGlobalIdTag, mHasMidNodesTag, mMaterialSetTag, mNeumannSetTag, mQaRecordTag, ncFile, NEUMANN_SET_TAG_NAME, moab::Interface::query_interface(), readMeshIface, reset(), and moab::Interface::tag_get_handle().

Referenced by factory().

◆ ~ReadNCDF()

moab::ReadNCDF::~ReadNCDF ( )
virtual

Destructor.

Definition at line 209 of file ReadNCDF.cpp.

210 {
212 }

References mdbImpl, readMeshIface, and moab::Interface::release_interface().

Member Function Documentation

◆ create_sideset_element()

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 1597 of file ReadNCDF.cpp.

1601 {
1602  // Get adjacent entities
1604  int to_dim = CN::Dimension( type );
1605  // for matching purposes , consider only corner nodes
1606  // basically, assume everything is matching if the corner nodes are matching, and
1607  // the number of total nodes is the same
1608  int nb_corner_nodes = CN::VerticesPerEntity( type );
1609  std::vector< EntityHandle > adj_ent;
1610  mdbImpl->get_adjacencies( &( connectivity[0] ), 1, to_dim, false, adj_ent );
1611 
1612  // For each entity, see if we can find a match
1613  // If we find a match, return it
1614  bool match_found = false;
1615  std::vector< EntityHandle > match_conn;
1616  // By default, assume positive sense
1617  sense = 1;
1618  for( unsigned int i = 0; i < adj_ent.size() && match_found == false; i++ )
1619  {
1620  // Get the connectivity
1621  error = mdbImpl->get_connectivity( &( adj_ent[i] ), 1, match_conn );
1622  if( error != MB_SUCCESS ) continue;
1623 
1624  // Make sure they have the same number of vertices (higher order elements ?)
1625  if( match_conn.size() != connectivity.size() ) continue;
1626 
1627  // Find a matching node
1628  std::vector< EntityHandle >::iterator iter;
1629  std::vector< EntityHandle >::iterator end_corner = match_conn.begin() + nb_corner_nodes;
1630  iter = std::find( match_conn.begin(), end_corner, connectivity[0] );
1631  if( iter == end_corner ) continue;
1632 
1633  // Rotate to match connectivity
1634  // rotate only corner nodes, ignore the rest from now on
1635  std::rotate( match_conn.begin(), iter, end_corner );
1636 
1637  bool they_match = true;
1638  for( int j = 1; j < nb_corner_nodes; j++ )
1639  {
1640  if( connectivity[j] != match_conn[j] )
1641  {
1642  they_match = false;
1643  break;
1644  }
1645  }
1646 
1647  // If we didn't get a match
1648  if( !they_match )
1649  {
1650  // Try the opposite sense
1651  they_match = true;
1652 
1653  int k = nb_corner_nodes - 1;
1654  for( int j = 1; j < nb_corner_nodes; )
1655  {
1656  if( connectivity[j] != match_conn[k] )
1657  {
1658  they_match = false;
1659  break;
1660  }
1661  ++j;
1662  --k;
1663  }
1664  // If they matched here, sense is reversed
1665  if( they_match ) sense = -1;
1666  }
1667  match_found = they_match;
1668  if( match_found ) handle = adj_ent[i];
1669  }
1670 
1671  // If we didn't find a match, create an element
1672  if( !match_found ) error = mdbImpl->create_element( type, &connectivity[0], connectivity.size(), handle );
1673 
1674  return error;
1675 }

References moab::Interface::create_element(), moab::CN::Dimension(), moab::error(), ErrorCode, moab::Interface::get_adjacencies(), moab::Interface::get_connectivity(), MB_SUCCESS, mdbImpl, and moab::CN::VerticesPerEntity().

Referenced by create_ss_elements().

◆ create_ss_elements()

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 1357 of file ReadNCDF.cpp.

1365 {
1366  // Determine entity type from element id
1367 
1368  // If there are dist. factors, create a vector to hold the array
1369  // and place this array as a tag onto the sideset meshset
1370 
1371  std::vector< double > temp_dist_factor_vector( num_dist_factors );
1372  std::vector< char > temp_string_storage( max_str_length + 1 );
1373  char* temp_string = &temp_string_storage[0];
1374  int temp_var;
1375  if( num_dist_factors )
1376  {
1377  INS_ID( temp_string, "dist_fact_ss%d", ss_seq_id );
1378  temp_var = -1;
1379  GET_1D_DBL_VAR( temp_string, temp_var, temp_dist_factor_vector );
1380  if( -1 == temp_var )
1381  {
1382  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting dist fact variable" );
1383  }
1384  }
1385 
1386  EntityType subtype;
1387  int num_side_nodes, num_elem_nodes;
1388  const EntityHandle* nodes;
1389  std::vector< EntityHandle > connectivity;
1390  int side_node_idx[32];
1391 
1392  int df_index = 0;
1393  int sense = 0;
1394  for( int i = 0; i < num_sides; i++ )
1395  {
1396  ExoIIElementType exoii_type;
1397  ReadBlockData block_data;
1398  block_data.elemType = EXOII_MAX_ELEM_TYPE;
1399 
1400  if( find_side_element_type( element_ids[i], exoii_type, block_data, df_index, side_list[i] ) != MB_SUCCESS )
1401  continue; // Isn't being read in this time
1402 
1403  EntityType type = ExoIIUtil::ExoIIElementMBEntity[exoii_type];
1404 
1405  EntityHandle ent_handle = element_ids[i] - block_data.startExoId + block_data.startMBId;
1406 
1407  const int side_num = side_list[i] - 1;
1408  if( type == MBHEX )
1409  {
1410  // Get the nodes of the element
1411  if( mdbImpl->get_connectivity( ent_handle, nodes, num_elem_nodes ) != MB_SUCCESS ) return MB_FAILURE;
1412 
1413  CN::SubEntityNodeIndices( type, num_elem_nodes, 2, side_num, subtype, num_side_nodes, side_node_idx );
1414  if( num_side_nodes <= 0 ) return MB_FAILURE;
1415 
1416  connectivity.resize( num_side_nodes );
1417  for( int k = 0; k < num_side_nodes; ++k )
1418  connectivity[k] = nodes[side_node_idx[k]];
1419 
1420  if( MB_SUCCESS != create_sideset_element( connectivity, subtype, ent_handle, sense ) ) return MB_FAILURE;
1421  if( 1 == sense )
1422  entities_to_add.push_back( ent_handle );
1423  else
1424  reverse_entities.push_back( ent_handle );
1425 
1426  // Read in distribution factor array
1427  if( num_dist_factors )
1428  {
1429  for( int k = 0; k < 4; k++ )
1430  dist_factor_vector.push_back( temp_dist_factor_vector[df_index++] );
1431  }
1432  }
1433 
1434  // If it is a Tet
1435  else if( type == MBTET )
1436  {
1437  // Get the nodes of the element
1438  if( mdbImpl->get_connectivity( ent_handle, nodes, num_elem_nodes ) != MB_SUCCESS ) return MB_FAILURE;
1439 
1440  CN::SubEntityNodeIndices( type, num_elem_nodes, 2, side_num, subtype, num_side_nodes, side_node_idx );
1441  if( num_side_nodes <= 0 ) return MB_FAILURE;
1442 
1443  connectivity.resize( num_side_nodes );
1444  for( int k = 0; k < num_side_nodes; ++k )
1445  connectivity[k] = nodes[side_node_idx[k]];
1446 
1447  if( MB_SUCCESS != create_sideset_element( connectivity, subtype, ent_handle, sense ) ) return MB_FAILURE;
1448  if( 1 == sense )
1449  entities_to_add.push_back( ent_handle );
1450  else
1451  reverse_entities.push_back( ent_handle );
1452 
1453  // Read in distribution factor array
1454  if( num_dist_factors )
1455  {
1456  for( int k = 0; k < 3; k++ )
1457  dist_factor_vector.push_back( temp_dist_factor_vector[df_index++] );
1458  }
1459  }
1460  else if( type == MBQUAD && exoii_type >= EXOII_SHELL && exoii_type <= EXOII_SHELL9 )
1461  {
1462  // ent_handle = CREATE_HANDLE(MBQUAD, base_id, error);
1463 
1464  // Just use this quad
1465  if( side_list[i] == 1 )
1466  {
1467  if( 1 == sense )
1468  entities_to_add.push_back( ent_handle );
1469  else
1470  reverse_entities.push_back( ent_handle );
1471 
1472  if( num_dist_factors )
1473  {
1474  for( int k = 0; k < 4; k++ )
1475  dist_factor_vector.push_back( temp_dist_factor_vector[df_index++] );
1476  }
1477 
1478  continue;
1479  }
1480  else if( side_list[i] == 2 )
1481  {
1482  reverse_entities.push_back( ent_handle );
1483 
1484  if( num_dist_factors )
1485  {
1486  for( int k = 0; k < 4; k++ )
1487  dist_factor_vector.push_back( temp_dist_factor_vector[df_index++] );
1488  }
1489 
1490  continue;
1491  }
1492  else
1493  {
1494  // Get the nodes of the element
1495  if( mdbImpl->get_connectivity( ent_handle, nodes, num_elem_nodes ) != MB_SUCCESS ) return MB_FAILURE;
1496 
1497  CN::SubEntityNodeIndices( type, num_elem_nodes, 1, side_num - 2, subtype, num_side_nodes,
1498  side_node_idx );
1499  if( num_side_nodes <= 0 ) return MB_FAILURE;
1500 
1501  connectivity.resize( num_side_nodes );
1502  for( int k = 0; k < num_side_nodes; ++k )
1503  connectivity[k] = nodes[side_node_idx[k]];
1504 
1505  if( MB_SUCCESS != create_sideset_element( connectivity, subtype, ent_handle, sense ) )
1506  return MB_FAILURE;
1507  if( 1 == sense )
1508  entities_to_add.push_back( ent_handle );
1509  else
1510  reverse_entities.push_back( ent_handle );
1511 
1512  if( num_dist_factors )
1513  {
1514  for( int k = 0; k < 2; k++ )
1515  dist_factor_vector.push_back( temp_dist_factor_vector[df_index++] );
1516  }
1517  }
1518  }
1519  // If it is a Quad
1520  else if( type == MBQUAD )
1521  {
1522  // Get the nodes of the element
1523  if( mdbImpl->get_connectivity( ent_handle, nodes, num_elem_nodes ) != MB_SUCCESS ) return MB_FAILURE;
1524 
1525  CN::SubEntityNodeIndices( type, num_elem_nodes, 1, side_num, subtype, num_side_nodes, side_node_idx );
1526  if( num_side_nodes <= 0 ) return MB_FAILURE;
1527 
1528  connectivity.resize( num_side_nodes );
1529  for( int k = 0; k < num_side_nodes; ++k )
1530  connectivity[k] = nodes[side_node_idx[k]];
1531 
1532  if( MB_SUCCESS != create_sideset_element( connectivity, subtype, ent_handle, sense ) ) return MB_FAILURE;
1533  if( 1 == sense )
1534  entities_to_add.push_back( ent_handle );
1535  else
1536  reverse_entities.push_back( ent_handle );
1537 
1538  // Read in distribution factor array
1539  if( num_dist_factors )
1540  {
1541  for( int k = 0; k < 2; k++ )
1542  dist_factor_vector.push_back( temp_dist_factor_vector[df_index++] );
1543  }
1544  }
1545  else if( type == MBTRI )
1546  {
1547  int side_offset = 0;
1548  if( number_dimensions() == 3 && side_list[i] <= 2 )
1549  {
1550  if( 1 == sense )
1551  entities_to_add.push_back( ent_handle );
1552  else
1553  reverse_entities.push_back( ent_handle );
1554  if( num_dist_factors )
1555  {
1556  for( int k = 0; k < 3; k++ )
1557  dist_factor_vector.push_back( temp_dist_factor_vector[df_index++] );
1558  }
1559  }
1560  else
1561  {
1562  if( number_dimensions() == 3 )
1563  {
1564  if( side_list[i] > 2 ) side_offset = 2;
1565  }
1566 
1567  // Get the nodes of the element
1568  if( mdbImpl->get_connectivity( ent_handle, nodes, num_elem_nodes ) != MB_SUCCESS ) return MB_FAILURE;
1569 
1570  CN::SubEntityNodeIndices( type, num_elem_nodes, 1, side_num - side_offset, subtype, num_side_nodes,
1571  side_node_idx );
1572  if( num_side_nodes <= 0 ) return MB_FAILURE;
1573 
1574  connectivity.resize( num_side_nodes );
1575  for( int k = 0; k < num_side_nodes; ++k )
1576  connectivity[k] = nodes[side_node_idx[k]];
1577 
1578  if( MB_SUCCESS != create_sideset_element( connectivity, subtype, ent_handle, sense ) )
1579  return MB_FAILURE;
1580  if( 1 == sense )
1581  entities_to_add.push_back( ent_handle );
1582  else
1583  reverse_entities.push_back( ent_handle );
1584 
1585  if( num_dist_factors )
1586  {
1587  for( int k = 0; k < 2; k++ )
1588  dist_factor_vector.push_back( temp_dist_factor_vector[df_index++] );
1589  }
1590  }
1591  }
1592  }
1593 
1594  return MB_SUCCESS;
1595 }

References create_sideset_element(), moab::ReadBlockData::elemType, moab::EXOII_MAX_ELEM_TYPE, moab::EXOII_SHELL, moab::EXOII_SHELL9, moab::ExoIIUtil::ExoIIElementMBEntity, find_side_element_type(), GET_1D_DBL_VAR, moab::Interface::get_connectivity(), INS_ID, max_str_length, MB_SET_ERR, MB_SUCCESS, MBHEX, MBQUAD, MBTET, MBTRI, mdbImpl, number_dimensions(), moab::ReadBlockData::startExoId, moab::ReadBlockData::startMBId, and moab::CN::SubEntityNodeIndices().

Referenced by read_sidesets().

◆ dimension_exists()

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

◆ exodus_file()

int moab::ReadNCDF::exodus_file ( )
private

exodus file bound to this object

◆ factory()

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

Definition at line 135 of file ReadNCDF.cpp.

136 {
137  return new ReadNCDF( iface );
138 }

References iface, and ReadNCDF().

Referenced by moab::ReaderWriterSet::ReaderWriterSet().

◆ find_side_element_type()

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 1677 of file ReadNCDF.cpp.

1682 {
1683  std::vector< ReadBlockData >::iterator iter, end_iter;
1684  iter = blocksLoading.begin();
1685  end_iter = blocksLoading.end();
1686  elem_type = EXOII_MAX_ELEM_TYPE;
1687 
1688  for( ; iter != end_iter; ++iter )
1689  {
1690  if( exodus_id >= ( *iter ).startExoId && exodus_id < ( *iter ).startExoId + ( *iter ).numElements )
1691  {
1692  elem_type = ( *iter ).elemType;
1693 
1694  // If we're not reading this block in
1695  if( !( *iter ).reading_in )
1696  {
1697  // Offset df_indes according to type
1698  if( elem_type >= EXOII_HEX && elem_type <= EXOII_HEX27 )
1699  df_index += 4;
1700  else if( elem_type >= EXOII_TETRA && elem_type <= EXOII_TETRA14 )
1701  df_index += 3;
1702  else if( elem_type >= EXOII_QUAD && elem_type <= EXOII_QUAD9 )
1703  df_index += 2;
1704  else if( elem_type >= EXOII_SHELL && elem_type <= EXOII_SHELL9 )
1705  {
1706  if( side_id == 1 || side_id == 2 )
1707  df_index += 4;
1708  else
1709  df_index += 2;
1710  }
1711  else if( elem_type >= EXOII_TRI && elem_type <= EXOII_TRI7 )
1712  df_index += 3;
1713 
1714  return MB_FAILURE;
1715  }
1716 
1717  block_data = *iter;
1718 
1719  return MB_SUCCESS;
1720  }
1721  }
1722  return MB_FAILURE;
1723 }

References blocksLoading, moab::EXOII_HEX, moab::EXOII_HEX27, moab::EXOII_MAX_ELEM_TYPE, moab::EXOII_QUAD, moab::EXOII_QUAD9, moab::EXOII_SHELL, moab::EXOII_SHELL9, moab::EXOII_TETRA, moab::EXOII_TETRA14, moab::EXOII_TRI, moab::EXOII_TRI7, and MB_SUCCESS.

Referenced by create_ss_elements().

◆ get_number_nodes()

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

◆ get_type() [1/2]

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

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

◆ get_type() [2/2]

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

◆ load_file()

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 284 of file ReadNCDF.cpp.

289 {
290  ErrorCode status;
291  int fail;
292 
293  int num_blocks = 0;
294  const int* blocks_to_load = 0;
295  if( subset_list )
296  {
297  if( subset_list->tag_list_length > 1 || !strcmp( subset_list->tag_list[0].tag_name, MATERIAL_SET_TAG_NAME ) )
298  {
299  MB_SET_ERR( MB_UNSUPPORTED_OPERATION, "ExodusII reader supports subset read only by material ID" );
300  }
301  if( subset_list->num_parts )
302  {
303  MB_SET_ERR( MB_UNSUPPORTED_OPERATION, "ExodusII reader does not support mesh partitioning" );
304  }
305  blocks_to_load = subset_list->tag_list[0].tag_values;
306  num_blocks = subset_list->tag_list[0].num_tag_values;
307  }
308 
309  // This function directs the reading of an exoii file, but doesn't do any of
310  // the actual work
311 
312  // See if opts has tdata.
313  ErrorCode rval;
314  std::string s;
315  rval = opts.get_str_option( "tdata", s );
316  if( MB_SUCCESS == rval && !s.empty() )
317  return update( exodus_file_name, opts, num_blocks, blocks_to_load, *file_set );
318 
319  reset();
320 
321  // 0. Open the file.
322 
323  // open netcdf/exodus file
324  fail = nc_open( exodus_file_name, 0, &ncFile );
325  if( NC_NOERR != fail )
326  {
327  MB_SET_ERR( MB_FILE_DOES_NOT_EXIST, "ReadNCDF:: problem opening Netcdf/Exodus II file " << exodus_file_name );
328  }
329 
330  // 1. Read the header
331  status = read_exodus_header();
332  if( MB_FAILURE == status ) return status;
333 
334  status = mdbImpl->get_entities_by_handle( 0, initRange );
335  if( MB_FAILURE == status ) return status;
336 
337  // 2. Read the nodes unless they've already been read before
338  status = read_nodes( file_id_tag );
339  if( MB_FAILURE == status ) return status;
340 
341  // 3.
342  // extra for polyhedra blocks
343  if( numberFaceBlocks_loading > 0 )
344  {
345  status = read_face_blocks_headers();
346  if( MB_FAILURE == status ) return status;
347  }
348 
349  status = read_block_headers( blocks_to_load, num_blocks );
350  if( MB_FAILURE == status ) return status;
351 
352  // 4. Read elements (might not read them, depending on active blocks)
353  if( numberFaceBlocks_loading > 0 )
354  {
355  status = read_polyhedra_faces();
356  if( MB_FAILURE == status ) return status;
357  }
358  status = read_elements( file_id_tag );
359  if( MB_FAILURE == status ) return status;
360 
361  // 5. Read global ids
362  status = read_global_ids();
363  if( MB_FAILURE == status ) return status;
364 
365  // 6. Read nodesets
366  status = read_nodesets();
367  if( MB_FAILURE == status ) return status;
368 
369  // 7. Read sidesets
370  status = read_sidesets();
371  if( MB_FAILURE == status ) return status;
372 
373  // 8. Read qa records
374  if( file_set )
375  {
376  status = read_qa_records( *file_set );
377  if( MB_FAILURE == status ) return status;
378  }
379 
380  // What about properties???
381 
382  // Close the file
383  fail = nc_close( ncFile );
384  if( NC_NOERR != fail )
385  {
386  MB_SET_ERR( MB_FAILURE, "Trouble closing file" );
387  }
388 
389  ncFile = 0;
390  return MB_SUCCESS;
391 }

References ErrorCode, moab::fail(), moab::Interface::get_entities_by_handle(), moab::FileOptions::get_str_option(), initRange, MATERIAL_SET_TAG_NAME, MB_FILE_DOES_NOT_EXIST, MB_SET_ERR, MB_SUCCESS, MB_UNSUPPORTED_OPERATION, mdbImpl, ncFile, moab::ReaderIface::SubsetList::num_parts, moab::ReaderIface::IDTag::num_tag_values, numberFaceBlocks_loading, read_block_headers(), read_elements(), read_exodus_header(), read_face_blocks_headers(), read_global_ids(), read_nodes(), read_nodesets(), read_polyhedra_faces(), read_qa_records(), read_sidesets(), reset(), moab::ReaderIface::SubsetList::tag_list, moab::ReaderIface::SubsetList::tag_list_length, moab::ReaderIface::IDTag::tag_name, moab::ReaderIface::IDTag::tag_values, and update().

◆ number_dimensions()

int moab::ReadNCDF::number_dimensions ( )
inlineprivate

number of dimensions in this exo file

Definition at line 255 of file ReadNCDF.hpp.

256 {
258 }

References numberDimensions_loading.

Referenced by create_ss_elements().

◆ read_block_headers()

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 553 of file ReadNCDF.cpp.

554 {
555  // Get the element block ids; keep this in a separate list,
556  // which is not offset by blockIdOffset; this list used later for
557  // reading block connectivity
558 
559  // Get the ids of all the blocks of this file we're reading in
560  std::vector< int > block_ids( numberElementBlocks_loading );
561  int nc_block_ids = -1;
562  GET_1D_INT_VAR( "eb_prop1", nc_block_ids, block_ids );
563  if( -1 == nc_block_ids )
564  {
565  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting eb_prop1 variable" );
566  }
567 
568  int exodus_id = 1;
569 
570  // If the active_block_id_list is NULL all blocks are active.
571  if( NULL == blocks_to_load || 0 == num_blocks )
572  {
573  blocks_to_load = &block_ids[0];
574  }
575 
576  std::vector< int > new_blocks( blocks_to_load, blocks_to_load + numberElementBlocks_loading );
577 
578  std::vector< int >::iterator iter, end_iter;
579  iter = block_ids.begin();
580  end_iter = block_ids.end();
581 
582  // Read header information and initialize header-type block information
583  int temp_dim;
584  std::vector< char > temp_string_storage( max_str_length + 1 );
585  char* temp_string = &temp_string_storage[0];
586  int block_seq_id = 1;
587 
588  for( ; iter != end_iter; ++iter, block_seq_id++ )
589  {
590  int num_elements;
591 
592  GET_DIMB( temp_dim, temp_string, "num_el_in_blk%d", block_seq_id, num_elements );
593 
594  // Don't read element type string for now, since it's an attrib
595  // on the connectivity
596  // Get the entity type corresponding to this ExoII element type
597  // ExoIIElementType elem_type =
598  // ExoIIUtil::static_element_name_to_type(element_type_string);
599 
600  // Tag each element block(mesh set) with enum for ElementType (SHELL, QUAD4, ....etc)
601  ReadBlockData block_data;
602  block_data.elemType = EXOII_MAX_ELEM_TYPE;
603  block_data.blockId = *iter;
604  block_data.startExoId = exodus_id;
605  block_data.numElements = num_elements;
606 
607  // If block is in 'blocks_to_load'----load it!
608  if( std::find( new_blocks.begin(), new_blocks.end(), *iter ) != new_blocks.end() )
609  block_data.reading_in = true;
610  else
611  block_data.reading_in = false;
612 
613  blocksLoading.push_back( block_data );
614  exodus_id += num_elements;
615  }
616 
617  return MB_SUCCESS;
618 }

References moab::ReadBlockData::blockId, blocksLoading, moab::ReadBlockData::elemType, moab::EXOII_MAX_ELEM_TYPE, GET_1D_INT_VAR, GET_DIMB, max_str_length, MB_SET_ERR, MB_SUCCESS, numberElementBlocks_loading, moab::ReadBlockData::numElements, moab::ReadBlockData::reading_in, and moab::ReadBlockData::startExoId.

Referenced by load_file(), and update().

◆ read_elements()

ErrorCode moab::ReadNCDF::read_elements ( const Tag file_id_tag)
private

read the element blocks

Definition at line 729 of file ReadNCDF.cpp.

730 {
731  // Read in elements
732 
733  int result = 0;
734 
735  // Initialize the nodeInLoadedBlocks vector
736  if( nodesInLoadedBlocks.size() < (size_t)( numberNodes_loading + 1 ) )
737  nodesInLoadedBlocks.resize( numberNodes_loading + 1 ); // this could be repeated?
738  memset( &nodesInLoadedBlocks[0], 0, ( numberNodes_loading + 1 ) * sizeof( char ) );
739 
740  std::vector< ReadBlockData >::iterator this_it;
741  this_it = blocksLoading.begin();
742 
743  int temp_dim;
744  std::vector< char > temp_string_storage( max_str_length + 1 );
745  char* temp_string = &temp_string_storage[0];
746 
747  int nc_var;
748  int block_seq_id = 1;
749  std::vector< int > dims;
750  size_t start[2] = { 0, 0 }, count[2];
751 
752  for( ; this_it != blocksLoading.end(); ++this_it, block_seq_id++ )
753  {
754  // If this block isn't to be read in --- continue
755  if( !( *this_it ).reading_in ) continue;
756 
757  // Get some information about this block
758  int block_id = ( *this_it ).blockId;
759  EntityHandle* conn = 0;
760 
761  // Get the ncdf connect variable and the element type
762  INS_ID( temp_string, "connect%d", block_seq_id );
763  GET_VAR( temp_string, nc_var, dims );
764  if( -1 == nc_var || 0 == nc_var )
765  { // try other var, for polyhedra blocks
766  // it could be polyhedra block, defined by fbconn and NFACED attribute
767  INS_ID( temp_string, "facconn%d", block_seq_id );
768  GET_VAR( temp_string, nc_var, dims );
769 
770  if( -1 == nc_var || 0 == nc_var )
771  {
772  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting connec or faccon variable" );
773  }
774  }
775  nc_type att_type;
776  size_t att_len;
777  int fail = nc_inq_att( ncFile, nc_var, "elem_type", &att_type, &att_len );
778  if( NC_NOERR != fail )
779  {
780  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting elem type attribute" );
781  }
782  std::vector< char > dum_str( att_len + 1 );
783  fail = nc_get_att_text( ncFile, nc_var, "elem_type", &dum_str[0] );
784  if( NC_NOERR != fail )
785  {
786  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting elem type" );
787  }
788  ExoIIElementType elem_type = ExoIIUtil::static_element_name_to_type( &dum_str[0] );
789  ( *this_it ).elemType = elem_type;
790 
791  int verts_per_element = ExoIIUtil::VerticesPerElement[( *this_it ).elemType];
792  int number_nodes = ( *this_it ).numElements * verts_per_element;
793  const EntityType mb_type = ExoIIUtil::ExoIIElementMBEntity[( *this_it ).elemType];
794 
795  if( mb_type == MBPOLYGON )
796  {
797  // need to read another variable, num_nod_per_el, and
798  int num_nodes_poly_conn;
799  GET_DIMB( temp_dim, temp_string, "num_nod_per_el%d", block_seq_id, num_nodes_poly_conn );
800  // read the connec1 array, which is of dimensions num_nodes_poly_conn (only one )
801  std::vector< int > connec;
802  connec.resize( num_nodes_poly_conn );
803  count[0] = num_nodes_poly_conn;
804 
805  fail = nc_get_vara_int( ncFile, nc_var, start, count, &connec[0] );
806  if( NC_NOERR != fail )
807  {
808  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting polygon connectivity " );
809  }
810  // ebepecnt1
811  std::vector< int > ebec;
812  ebec.resize( this_it->numElements );
813  // Get the ncdf connect variable and the element type
814  INS_ID( temp_string, "ebepecnt%d", block_seq_id );
815  GET_VAR( temp_string, nc_var, dims );
816  count[0] = this_it->numElements;
817  fail = nc_get_vara_int( ncFile, nc_var, start, count, &ebec[0] );
818  if( NC_NOERR != fail )
819  {
820  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting polygon nodes per elem " );
821  }
822  EntityHandle ms_handle;
824  return MB_FAILURE;
825  // create polygons one by one, and put them in the list
826  // also need to read the number of nodes for each polygon, then create
827  std::vector< EntityHandle > polyConn;
828  int ix = 0;
829  for( int i = 0; i < this_it->numElements; i++ )
830  {
831  polyConn.resize( ebec[i] );
832  for( int j = 0; j < ebec[i]; j++ )
833  {
834  nodesInLoadedBlocks[connec[ix]] = 1;
835  polyConn[j] = vertexOffset + connec[ix++];
836  }
837  EntityHandle newp;
838  /*
839  * ErrorCode create_element(const EntityType type,
840  const EntityHandle *connectivity,
841  const int num_vertices,
842  EntityHandle &element_handle)
843  */
844  if( mdbImpl->create_element( MBPOLYGON, &polyConn[0], ebec[i], newp ) != MB_SUCCESS ) return MB_FAILURE;
845 
846  // add the element in order
847  this_it->polys.push_back( newp );
848  if( mdbImpl->add_entities( ms_handle, &newp, 1 ) != MB_SUCCESS ) return MB_FAILURE;
849  }
850  // Set the block id with an offset
851  if( mdbImpl->tag_set_data( mMaterialSetTag, &ms_handle, 1, &block_id ) != MB_SUCCESS ) return MB_FAILURE;
852  if( mdbImpl->tag_set_data( mGlobalIdTag, &ms_handle, 1, &block_id ) != MB_SUCCESS ) return MB_FAILURE;
853  }
854  else if( mb_type == MBPOLYHEDRON )
855  {
856  // need to read another variable, num_fac_per_el
857  int num_fac_per_el;
858  GET_DIMB( temp_dim, temp_string, "num_fac_per_el%d", block_seq_id, num_fac_per_el );
859  // read the fbconn array, which is of dimensions num_nod_per_fa (only one dimension)
860  std::vector< int > facconn;
861  facconn.resize( num_fac_per_el );
862  count[0] = num_fac_per_el;
863 
864  fail = nc_get_vara_int( ncFile, nc_var, start, count, &facconn[0] );
865  if( NC_NOERR != fail )
866  {
867  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting polygon connectivity " );
868  }
869  // ebepecnt1
870  std::vector< int > ebepecnt;
871  ebepecnt.resize( this_it->numElements );
872  // Get the ncdf connect variable and the element type
873  INS_ID( temp_string, "ebepecnt%d", block_seq_id );
874  GET_VAR( temp_string, nc_var, dims );
875  count[0] = this_it->numElements;
876  fail = nc_get_vara_int( ncFile, nc_var, start, count, &ebepecnt[0] );
877  if( NC_NOERR != fail )
878  {
879  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting polygon nodes per elem " );
880  }
881  EntityHandle ms_handle;
883  return MB_FAILURE;
884  // create polygons one by one, and put them in the list
885  // also need to read the number of nodes for each polygon, then create
886  std::vector< EntityHandle > polyConn;
887  int ix = 0;
888  for( int i = 0; i < this_it->numElements; i++ )
889  {
890  polyConn.resize( ebepecnt[i] );
891  for( int j = 0; j < ebepecnt[i]; j++ )
892  {
893  polyConn[j] = polyfaces[facconn[ix++] - 1];
894  }
895  EntityHandle newp;
896  /*
897  * ErrorCode create_element(const EntityType type,
898  const EntityHandle *connectivity,
899  const int num_vertices,
900  EntityHandle &element_handle)
901  */
902  if( mdbImpl->create_element( MBPOLYHEDRON, &polyConn[0], ebepecnt[i], newp ) != MB_SUCCESS )
903  return MB_FAILURE;
904 
905  // add the element in order
906  this_it->polys.push_back( newp );
907  if( mdbImpl->add_entities( ms_handle, &newp, 1 ) != MB_SUCCESS ) return MB_FAILURE;
908  }
909  // Set the block id with an offset
910  if( mdbImpl->tag_set_data( mMaterialSetTag, &ms_handle, 1, &block_id ) != MB_SUCCESS ) return MB_FAILURE;
911  if( mdbImpl->tag_set_data( mGlobalIdTag, &ms_handle, 1, &block_id ) != MB_SUCCESS ) return MB_FAILURE;
912  }
913  else // this is regular
914  {
915  // Allocate an array to read in connectivity data
916  readMeshIface->get_element_connect( this_it->numElements, verts_per_element, mb_type, this_it->startExoId,
917  this_it->startMBId, conn );
918 
919  // Create a range for this sequence of elements
920  EntityHandle start_range, end_range;
921  start_range = ( *this_it ).startMBId;
922  end_range = start_range + ( *this_it ).numElements - 1;
923 
924  Range new_range( start_range, end_range );
925  // Range<EntityHandle> new_range((*this_it).startMBId,
926  // (*this_it).startMBId + (*this_it).numElements - 1);
927 
928  // Create a MBSet for this block and set the material tag
929 
930  EntityHandle ms_handle;
932  return MB_FAILURE;
933 
934  if( mdbImpl->add_entities( ms_handle, new_range ) != MB_SUCCESS ) return MB_FAILURE;
935 
936  int mid_nodes[4];
937  CN::HasMidNodes( mb_type, verts_per_element, mid_nodes );
938  if( mdbImpl->tag_set_data( mHasMidNodesTag, &ms_handle, 1, mid_nodes ) != MB_SUCCESS ) return MB_FAILURE;
939 
940  // Just a check because the following code won't work if this case fails
941  assert( sizeof( EntityHandle ) >= sizeof( int ) );
942 
943  // tmp_ptr is of type int* and points at the same place as conn
944  int* tmp_ptr = reinterpret_cast< int* >( conn );
945 
946  // Read the connectivity into that memory, this will take only part of the array
947  // 1/2 if sizeof(EntityHandle) == 64 bits.
948  count[0] = this_it->numElements;
949  count[1] = verts_per_element;
950  fail = nc_get_vara_int( ncFile, nc_var, start, count, tmp_ptr );
951  if( NC_NOERR != fail )
952  {
953  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting connectivity" );
954  }
955 
956  // Convert from exodus indices to vertex handles.
957  // Iterate backwards in case handles are larger than ints.
958  for( int i = number_nodes - 1; i >= 0; --i )
959  {
960  if( (unsigned)tmp_ptr[i] >= nodesInLoadedBlocks.size() )
961  {
962  MB_SET_ERR( MB_FAILURE, "Invalid node ID in block connectivity" );
963  }
964  nodesInLoadedBlocks[tmp_ptr[i]] = 1;
965  conn[i] = static_cast< EntityHandle >( tmp_ptr[i] ) + vertexOffset;
966  }
967 
968  // Adjust connectivity order if necessary
969  const int* reorder = exodus_elem_order_map[mb_type][verts_per_element];
970  if( reorder ) ReadUtilIface::reorder( reorder, conn, this_it->numElements, verts_per_element );
971 
972  readMeshIface->update_adjacencies( ( *this_it ).startMBId, ( *this_it ).numElements,
973  ExoIIUtil::VerticesPerElement[( *this_it ).elemType], conn );
974 
975  if( result == -1 )
976  {
977  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: error getting element connectivity for block " << block_id );
978  }
979 
980  // Set the block id with an offset
981  if( mdbImpl->tag_set_data( mMaterialSetTag, &ms_handle, 1, &block_id ) != MB_SUCCESS ) return MB_FAILURE;
982  if( mdbImpl->tag_set_data( mGlobalIdTag, &ms_handle, 1, &block_id ) != MB_SUCCESS ) return MB_FAILURE;
983 
984  if( file_id_tag )
985  {
986  Range range;
987  range.insert( this_it->startMBId, this_it->startMBId + this_it->numElements - 1 );
988  readMeshIface->assign_ids( *file_id_tag, range, this_it->startExoId );
989  }
990  } // end regular block (not polygon)
991  }
992 
993  return MB_SUCCESS;
994 }

References moab::Interface::add_entities(), moab::ReadUtilIface::assign_ids(), blocksLoading, moab::Interface::create_element(), moab::Interface::create_meshset(), exodus_elem_order_map, moab::ExoIIUtil::ExoIIElementMBEntity, moab::fail(), GET_DIMB, moab::ReadUtilIface::get_element_connect(), GET_VAR, moab::CN::HasMidNodes(), INS_ID, moab::Range::insert(), max_str_length, MB_SET_ERR, MB_SUCCESS, MBPOLYGON, MBPOLYHEDRON, mdbImpl, MESHSET_SET, MESHSET_TRACK_OWNER, mGlobalIdTag, mHasMidNodesTag, mMaterialSetTag, ncFile, nodesInLoadedBlocks, numberNodes_loading, polyfaces, readMeshIface, moab::ReadUtilIface::reorder(), moab::ExoIIUtil::static_element_name_to_type(), moab::Interface::tag_set_data(), moab::ReadUtilIface::update_adjacencies(), vertexOffset, and moab::ExoIIUtil::VerticesPerElement.

Referenced by load_file().

◆ read_exodus_header()

ErrorCode moab::ReadNCDF::read_exodus_header ( )
private

read the header from the ExoII file

Definition at line 393 of file ReadNCDF.cpp.

394 {
395  CPU_WORD_SIZE = sizeof( double ); // With ExodusII version 2, all floats
396  IO_WORD_SIZE = sizeof( double ); // should be changed to doubles
397 
398  // NetCDF doesn't check its own limits on file read, so check
399  // them here so it doesn't corrupt memory any more than absolutely
400  // necessary.
401  int ndims;
402  int fail = nc_inq_ndims( ncFile, &ndims );
403  if( NC_NOERR != fail || ndims > NC_MAX_DIMS )
404  {
405  MB_SET_ERR( MB_FAILURE, "ReadNCDF: File contains " << ndims << " dims but NetCDF library supports only "
406  << (int)NC_MAX_DIMS );
407  }
408  int nvars;
409  fail = nc_inq_nvars( ncFile, &nvars );
410  if( nvars > NC_MAX_VARS )
411  {
412  MB_SET_ERR( MB_FAILURE, "ReadNCDF: File contains " << nvars << " vars but NetCDF library supports only "
413  << (int)NC_MAX_VARS );
414  }
415 
416  // Get the attributes
417 
418  // Get the word size, scalar value
419  nc_type att_type;
420  size_t att_len;
421  fail = nc_inq_att( ncFile, NC_GLOBAL, "floating_point_word_size", &att_type, &att_len );
422  if( NC_NOERR != fail )
423  {
424  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting floating_point_word_size attribute" );
425  }
426  if( att_type != NC_INT || att_len != 1 )
427  {
428  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Word size didn't have type int or size 1" );
429  }
430  fail = nc_get_att_int( ncFile, NC_GLOBAL, "floating_point_word_size", &IO_WORD_SIZE );
431  if( NC_NOERR != fail )
432  {
433  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Trouble getting word size" );
434  }
435 
436  // Exodus version
437  fail = nc_inq_att( ncFile, NC_GLOBAL, "version", &att_type, &att_len );
438  if( NC_NOERR != fail )
439  {
440  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting version attribute" );
441  }
442  if( att_type != NC_FLOAT || att_len != 1 )
443  {
444  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Version didn't have type float or size 1" );
445  }
446  // float version = temp_att->as_float(0);
447 
448  // Read in initial variables
449  int temp_dim;
450  GET_DIM( temp_dim, "num_dim", numberDimensions_loading );
451  GET_DIM( temp_dim, "num_nodes", numberNodes_loading );
452  GET_DIM( temp_dim, "num_elem", numberElements_loading );
453  GET_DIM( temp_dim, "num_el_blk", numberElementBlocks_loading );
454  GET_DIM( temp_dim, "num_fa_blk", numberFaceBlocks_loading ); // for polyhedra blocks
455  GET_DIM( temp_dim, "num_elem", numberElements_loading );
456  GET_DIM( temp_dim, "num_node_sets", numberNodeSets_loading );
457  GET_DIM( temp_dim, "num_side_sets", numberSideSets_loading );
458  GET_DIM( temp_dim, "len_string", max_str_length );
459  GET_DIM( temp_dim, "len_line", max_line_length );
460 
461  // Title; why are we even bothering if we're not going to keep it???
462  fail = nc_inq_att( ncFile, NC_GLOBAL, "title", &att_type, &att_len );
463  if( NC_NOERR != fail )
464  {
465  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting title attribute" );
466  }
467  if( att_type != NC_CHAR )
468  {
469  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: title didn't have type char" );
470  }
471  char* title = new char[att_len + 1];
472  fail = nc_get_att_text( ncFile, NC_GLOBAL, "title", title );
473  if( NC_NOERR != fail )
474  {
475  delete[] title;
476  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: trouble getting title" );
477  }
478  delete[] title;
479 
480  return MB_SUCCESS;
481 }

References CPU_WORD_SIZE, moab::fail(), GET_DIM, IO_WORD_SIZE, max_line_length, max_str_length, MB_SET_ERR, MB_SUCCESS, ncFile, numberDimensions_loading, numberElementBlocks_loading, numberElements_loading, numberFaceBlocks_loading, numberNodes_loading, numberNodeSets_loading, and numberSideSets_loading.

Referenced by load_file(), read_tag_values(), and update().

◆ read_face_blocks_headers()

ErrorCode moab::ReadNCDF::read_face_blocks_headers ( )
private

Definition at line 620 of file ReadNCDF.cpp.

621 {
622  // Get the ids of all the blocks of this file we're reading in
623  std::vector< int > fblock_ids( numberFaceBlocks_loading );
624  int nc_fblock_ids = -1;
625  GET_1D_INT_VAR( "fa_prop1", nc_fblock_ids, fblock_ids );
626  if( -1 == nc_fblock_ids )
627  {
628  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting fa_prop1 variable" );
629  }
630 
631  int temp_dim;
632  std::vector< char > temp_string_storage( max_str_length + 1 );
633  char* temp_string = &temp_string_storage[0];
634  // int block_seq_id = 1;
635  int exodus_id = 1;
636 
637  for( int i = 1; i <= numberFaceBlocks_loading; i++ )
638  {
639  int num_elements;
640 
641  GET_DIMB( temp_dim, temp_string, "num_fa_in_blk%d", i, num_elements );
642 
643  // Don't read element type string for now, since it's an attrib
644  // on the connectivity
645  // Get the entity type corresponding to this ExoII element type
646  // ExoIIElementType elem_type =
647  // ExoIIUtil::static_element_name_to_type(element_type_string);
648 
649  // Tag each element block(mesh set) with enum for ElementType (SHELL, QUAD4, ....etc)
650  ReadFaceBlockData fblock_data;
651  fblock_data.faceBlockId = i;
652  fblock_data.startExoId = exodus_id;
653  fblock_data.numElements = num_elements;
654 
655  faceBlocksLoading.push_back( fblock_data );
656  exodus_id += num_elements;
657  }
658  return MB_SUCCESS;
659 }

References moab::ReadFaceBlockData::faceBlockId, faceBlocksLoading, GET_1D_INT_VAR, GET_DIMB, max_str_length, MB_SET_ERR, MB_SUCCESS, numberFaceBlocks_loading, moab::ReadFaceBlockData::numElements, and moab::ReadFaceBlockData::startExoId.

Referenced by load_file().

◆ read_global_ids()

ErrorCode moab::ReadNCDF::read_global_ids ( )
private

read in the global element ids

Definition at line 996 of file ReadNCDF.cpp.

997 {
998  // Read in the map from the exodus file
999  std::vector< int > ids( std::max( numberElements_loading, numberNodes_loading ) );
1000 
1001  int varid = -1;
1002  GET_1D_INT_VAR( "elem_map", varid, ids );
1003  if( -1 != varid )
1004  {
1005  std::vector< ReadBlockData >::iterator iter;
1006  int id_pos = 0;
1007  for( iter = blocksLoading.begin(); iter != blocksLoading.end(); ++iter )
1008  {
1009  if( iter->reading_in )
1010  {
1011  if( iter->elemType != EXOII_POLYGON && iter->elemType != EXOII_POLYHEDRON )
1012  {
1013  if( iter->startMBId != 0 )
1014  {
1015  Range range( iter->startMBId, iter->startMBId + iter->numElements - 1 );
1016  ErrorCode error = mdbImpl->tag_set_data( mGlobalIdTag, range, &ids[id_pos] );
1017  if( error != MB_SUCCESS ) return error;
1018  id_pos += iter->numElements;
1019  }
1020  else
1021  return MB_FAILURE;
1022  }
1023  else // polygons or polyhedrons; use id from elements
1024  {
1025  for( std::vector< EntityHandle >::iterator eit = iter->polys.begin(); eit != iter->polys.end();
1026  eit++, id_pos++ )
1027  {
1028  EntityHandle peh = *eit;
1029  if( mdbImpl->tag_set_data( mGlobalIdTag, &peh, 1, &ids[id_pos] ) != MB_SUCCESS )
1030  return MB_FAILURE;
1031  }
1032  }
1033  }
1034  }
1035  }
1036 
1037  // Read in node map next
1038  varid = -1;
1039  GET_1D_INT_VAR( "node_num_map", varid, ids );
1040  if( -1 != varid )
1041  {
1043  ErrorCode error = mdbImpl->tag_set_data( mGlobalIdTag, range, &ids[0] );
1044  if( MB_SUCCESS != error ) MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem setting node global ids" );
1045  }
1046 
1047  return MB_SUCCESS;
1048 }

References blocksLoading, moab::error(), ErrorCode, moab::EXOII_POLYGON, moab::EXOII_POLYHEDRON, GET_1D_INT_VAR, MB_SET_ERR, MB_START_ID, MB_SUCCESS, mdbImpl, mGlobalIdTag, numberElements_loading, numberNodes_loading, moab::Interface::tag_set_data(), and vertexOffset.

Referenced by load_file().

◆ read_nodes()

ErrorCode moab::ReadNCDF::read_nodes ( const Tag file_id_tag)
private

read the nodes

Definition at line 483 of file ReadNCDF.cpp.

484 {
485  // Read the nodes into memory
486 
487  assert( 0 != ncFile );
488 
489  // Create a sequence to hold the node coordinates
490  // Get the current number of entities and start at the next slot
491 
492  EntityHandle node_handle = 0;
493  std::vector< double* > arrays;
494  readMeshIface->get_node_coords( 3, numberNodes_loading, MB_START_ID, node_handle, arrays );
495 
496  vertexOffset = ID_FROM_HANDLE( node_handle ) - MB_START_ID;
497 
498  // Read in the coordinates
499  int fail;
500  int coord = 0;
501  nc_inq_varid( ncFile, "coord", &coord );
502 
503  // Single var for all coords
504  size_t start[2] = { 0, 0 }, count[2] = { 1, static_cast< size_t >( numberNodes_loading ) };
505  if( coord )
506  {
507  for( int d = 0; d < numberDimensions_loading; ++d )
508  {
509  start[0] = d;
510  fail = nc_get_vara_double( ncFile, coord, start, count, arrays[d] );
511  if( NC_NOERR != fail )
512  {
513  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting " << (char)( 'x' + d ) << " coord array" );
514  }
515  }
516  }
517  // Var for each coord
518  else
519  {
520  char varname[] = "coord ";
521  for( int d = 0; d < numberDimensions_loading; ++d )
522  {
523  varname[5] = 'x' + (char)d;
524  fail = nc_inq_varid( ncFile, varname, &coord );
525  if( NC_NOERR != fail )
526  {
527  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting " << (char)( 'x' + d ) << " coord variable" );
528  }
529 
530  fail = nc_get_vara_double( ncFile, coord, start, &count[1], arrays[d] );
531  if( NC_NOERR != fail )
532  {
533  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting " << (char)( 'x' + d ) << " coord array" );
534  }
535  }
536  }
537 
538  // Zero out any coord values that are in database but not in file
539  // (e.g. if MOAB has 3D coords but file is 2D then set Z coords to zero.)
540  for( unsigned d = numberDimensions_loading; d < arrays.size(); ++d )
541  std::fill( arrays[d], arrays[d] + numberNodes_loading, 0.0 );
542 
543  if( file_id_tag )
544  {
545  Range nodes;
546  nodes.insert( node_handle, node_handle + numberNodes_loading - 1 );
547  readMeshIface->assign_ids( *file_id_tag, nodes, vertexOffset );
548  }
549 
550  return MB_SUCCESS;
551 }

References moab::ReadUtilIface::assign_ids(), moab::fail(), moab::ReadUtilIface::get_node_coords(), moab::ID_FROM_HANDLE(), moab::Range::insert(), MB_SET_ERR, MB_START_ID, MB_SUCCESS, ncFile, numberDimensions_loading, numberNodes_loading, readMeshIface, and vertexOffset.

Referenced by load_file().

◆ read_nodesets()

ErrorCode moab::ReadNCDF::read_nodesets ( )
private

read the nodesets into meshsets

Definition at line 1050 of file ReadNCDF.cpp.

1051 {
1052  // Read in the nodesets for the model
1053 
1054  if( 0 == numberNodeSets_loading ) return MB_SUCCESS;
1055  std::vector< int > id_array( numberNodeSets_loading );
1056 
1057  // Read in the nodeset ids
1058  int nc_var;
1059  GET_1D_INT_VAR( "ns_prop1", nc_var, id_array );
1060  if( -1 == nc_var )
1061  {
1062  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting ns_prop1 variable" );
1063  }
1064 
1065  // Use a vector of ints to read node handles
1066  std::vector< int > node_handles;
1067 
1068  int i;
1069  std::vector< char > temp_string_storage( max_str_length + 1 );
1070  char* temp_string = &temp_string_storage[0];
1071  int temp_dim;
1072  for( i = 0; i < numberNodeSets_loading; i++ )
1073  {
1074  // Get nodeset parameters
1075  int number_nodes_in_set;
1076  int number_dist_factors_in_set;
1077 
1078  GET_DIMB( temp_dim, temp_string, "num_nod_ns%d", i + 1, number_nodes_in_set );
1079  GET_DIMB( temp_dim, temp_string, "num_df_ns%d", i + 1, number_dist_factors_in_set );
1080 
1081  // Need to new a vector to store dist. factors
1082  // This vector * gets stored as a tag on the sideset meshset
1083  std::vector< double > temp_dist_factor_vector( number_nodes_in_set );
1084  if( number_dist_factors_in_set != 0 )
1085  {
1086  INS_ID( temp_string, "dist_fact_ns%d", i + 1 );
1087  GET_1D_DBL_VAR( temp_string, temp_dim, temp_dist_factor_vector );
1088  if( -1 == temp_dim )
1089  {
1090  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting dist fact variable" );
1091  }
1092  }
1093 
1094  // Size new arrays and get ids and distribution factors
1095  if( node_handles.size() < (unsigned int)number_nodes_in_set )
1096  {
1097  node_handles.reserve( number_nodes_in_set );
1098  node_handles.resize( number_nodes_in_set );
1099  }
1100 
1101  INS_ID( temp_string, "node_ns%d", i + 1 );
1102  int temp_var = -1;
1103  GET_1D_INT_VAR( temp_string, temp_var, node_handles );
1104  if( -1 == temp_var )
1105  {
1106  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting nodeset node variable" );
1107  }
1108 
1109  // Maybe there is already a nodesets meshset here we can append to
1110  Range child_meshsets;
1111  if( mdbImpl->get_entities_by_handle( 0, child_meshsets ) != MB_SUCCESS ) return MB_FAILURE;
1112 
1113  child_meshsets = subtract( child_meshsets, initRange );
1114 
1115  Range::iterator iter, end_iter;
1116  iter = child_meshsets.begin();
1117  end_iter = child_meshsets.end();
1118 
1119  EntityHandle ns_handle = 0;
1120  for( ; iter != end_iter; ++iter )
1121  {
1122  int nodeset_id;
1123  if( mdbImpl->tag_get_data( mDirichletSetTag, &( *iter ), 1, &nodeset_id ) != MB_SUCCESS ) continue;
1124 
1125  if( id_array[i] == nodeset_id )
1126  {
1127  // Found the meshset
1128  ns_handle = *iter;
1129  break;
1130  }
1131  }
1132 
1133  std::vector< EntityHandle > nodes_of_nodeset;
1134  if( ns_handle )
1135  if( mdbImpl->get_entities_by_handle( ns_handle, nodes_of_nodeset, true ) != MB_SUCCESS ) return MB_FAILURE;
1136 
1137  // Make these into entity handles
1138  // TODO: could we have read it into EntityHandle sized array in the first place?
1139  int j, temp;
1140  std::vector< EntityHandle > nodes;
1141  std::vector< double > dist_factor_vector;
1142  for( j = 0; j < number_nodes_in_set; j++ )
1143  {
1144  // See if this node is one we're currently reading in
1145  if( nodesInLoadedBlocks[node_handles[j]] == 1 )
1146  {
1147  // Make sure that it already isn't in a nodeset
1148  unsigned int node_id = CREATE_HANDLE( MBVERTEX, node_handles[j] + vertexOffset, temp );
1149  if( !ns_handle ||
1150  std::find( nodes_of_nodeset.begin(), nodes_of_nodeset.end(), node_id ) == nodes_of_nodeset.end() )
1151  {
1152  nodes.push_back( node_id );
1153 
1154  if( number_dist_factors_in_set != 0 ) dist_factor_vector.push_back( temp_dist_factor_vector[j] );
1155  }
1156  }
1157  }
1158 
1159  // No nodes to add
1160  if( nodes.empty() ) continue;
1161 
1162  // If there was no meshset found --> create one
1163  if( ns_handle == 0 )
1164  {
1165  if( mdbImpl->create_meshset( MESHSET_ORDERED | MESHSET_TRACK_OWNER, ns_handle ) != MB_SUCCESS )
1166  return MB_FAILURE;
1167 
1168  // Set a tag signifying dirichlet bc
1169  // TODO: create this tag another way
1170 
1171  int nodeset_id = id_array[i];
1172  if( mdbImpl->tag_set_data( mDirichletSetTag, &ns_handle, 1, &nodeset_id ) != MB_SUCCESS ) return MB_FAILURE;
1173  if( mdbImpl->tag_set_data( mGlobalIdTag, &ns_handle, 1, &nodeset_id ) != MB_SUCCESS ) return MB_FAILURE;
1174 
1175  if( !dist_factor_vector.empty() )
1176  {
1177  int size = dist_factor_vector.size();
1178  const void* data = &dist_factor_vector[0];
1179  if( mdbImpl->tag_set_by_ptr( mDistFactorTag, &ns_handle, 1, &data, &size ) != MB_SUCCESS )
1180  return MB_FAILURE;
1181  }
1182  }
1183  else if( !dist_factor_vector.empty() )
1184  {
1185  // Append dist factors to vector
1186  const void* ptr = 0;
1187  int size = 0;
1188  if( mdbImpl->tag_get_by_ptr( mDistFactorTag, &ns_handle, 1, &ptr, &size ) != MB_SUCCESS ) return MB_FAILURE;
1189 
1190  const double* data = reinterpret_cast< const double* >( ptr );
1191  dist_factor_vector.reserve( dist_factor_vector.size() + size );
1192  std::copy( data, data + size, std::back_inserter( dist_factor_vector ) );
1193  size = dist_factor_vector.size();
1194  ptr = &dist_factor_vector[0];
1195  if( mdbImpl->tag_set_by_ptr( mDistFactorTag, &ns_handle, 1, &ptr, &size ) != MB_SUCCESS ) return MB_FAILURE;
1196  }
1197 
1198  // Add the nodes to the meshset
1199  if( mdbImpl->add_entities( ns_handle, &nodes[0], nodes.size() ) != MB_SUCCESS ) return MB_FAILURE;
1200  }
1201 
1202  return MB_SUCCESS;
1203 }

References moab::Interface::add_entities(), moab::Range::begin(), moab::CREATE_HANDLE(), moab::Interface::create_meshset(), moab::Range::end(), GET_1D_DBL_VAR, GET_1D_INT_VAR, GET_DIMB, moab::Interface::get_entities_by_handle(), initRange, INS_ID, max_str_length, MB_SET_ERR, MB_SUCCESS, MBVERTEX, mdbImpl, mDirichletSetTag, mDistFactorTag, MESHSET_TRACK_OWNER, mGlobalIdTag, nodesInLoadedBlocks, numberNodeSets_loading, size, moab::subtract(), moab::Interface::tag_get_by_ptr(), moab::Interface::tag_get_data(), moab::Interface::tag_set_by_ptr(), moab::Interface::tag_set_data(), and vertexOffset.

Referenced by load_file().

◆ read_polyhedra_faces()

ErrorCode moab::ReadNCDF::read_polyhedra_faces ( )
private

Definition at line 660 of file ReadNCDF.cpp.

661 {
662  int temp_dim;
663  std::vector< char > temp_string_storage( max_str_length + 1 );
664  char* temp_string = &temp_string_storage[0];
666  size_t start[1] = { 0 }, count[1]; // one dim arrays only here!
667  std::vector< ReadFaceBlockData >::iterator this_it;
668  this_it = faceBlocksLoading.begin();
669 
670  int fblock_seq_id = 1;
671  std::vector< int > dims;
672  int nc_var, fail;
673 
674  for( ; this_it != faceBlocksLoading.end(); ++this_it, fblock_seq_id++ )
675  {
676  int num_fa_in_blk;
677  GET_DIMB( temp_dim, temp_string, "num_fa_in_blk%d", fblock_seq_id, num_fa_in_blk );
678  int num_nod_per_fa;
679  GET_DIMB( temp_dim, temp_string, "num_nod_per_fa%d", fblock_seq_id, num_nod_per_fa );
680  // Get the ncdf connect variable and the element type
681  INS_ID( temp_string, "fbconn%d", fblock_seq_id );
682  GET_VAR( temp_string, nc_var, dims );
683  std::vector< int > fbconn;
684  fbconn.resize( num_nod_per_fa );
685  count[0] = num_nod_per_fa;
686 
687  fail = nc_get_vara_int( ncFile, nc_var, start, count, &fbconn[0] );
688  if( NC_NOERR != fail )
689  {
690  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting face polyhedra connectivity " );
691  }
692  std::vector< int > fbepecnt;
693  fbepecnt.resize( num_fa_in_blk );
694  INS_ID( temp_string, "fbepecnt%d", fblock_seq_id );
695  GET_VAR( temp_string, nc_var, dims );
696  count[0] = num_fa_in_blk;
697  fail = nc_get_vara_int( ncFile, nc_var, start, count, &fbepecnt[0] );
698  if( NC_NOERR != fail )
699  {
700  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting face polyhedra connectivity " );
701  }
702  // now we can create some polygons
703  std::vector< EntityHandle > polyConn;
704  int ix = 0;
705  for( int i = 0; i < num_fa_in_blk; i++ )
706  {
707  polyConn.resize( fbepecnt[i] );
708  for( int j = 0; j < fbepecnt[i]; j++ )
709  {
710  nodesInLoadedBlocks[fbconn[ix]] = 1;
711  polyConn[j] = vertexOffset + fbconn[ix++];
712  }
713  EntityHandle newp;
714  /*
715  * ErrorCode create_element(const EntityType type,
716  const EntityHandle *connectivity,
717  const int num_vertices,
718  EntityHandle &element_handle)
719  */
720  if( mdbImpl->create_element( MBPOLYGON, &polyConn[0], fbepecnt[i], newp ) != MB_SUCCESS ) return MB_FAILURE;
721 
722  // add the element in order
723  polyfaces.push_back( newp );
724  }
725  }
726  return MB_SUCCESS;
727 }

References moab::Interface::create_element(), faceBlocksLoading, moab::fail(), GET_DIMB, GET_VAR, INS_ID, max_str_length, MB_SET_ERR, MB_SUCCESS, MBPOLYGON, mdbImpl, ncFile, nodesInLoadedBlocks, numberNodes_loading, polyfaces, and vertexOffset.

Referenced by load_file().

◆ read_qa_information()

ErrorCode moab::ReadNCDF::read_qa_information ( std::vector< std::string > &  qa_record_list)
private

Definition at line 1751 of file ReadNCDF.cpp.

1752 {
1753  // Inquire on the genesis file to find the number of qa records
1754 
1755  int number_records = 0;
1756 
1757  int temp_dim;
1758  GET_DIM( temp_dim, "num_qa_rec", number_records );
1759  std::vector< char > data( max_str_length + 1 );
1760 
1761  for( int i = 0; i < number_records; i++ )
1762  {
1763  for( int j = 0; j < 4; j++ )
1764  {
1765  data[max_str_length] = '\0';
1766  if( read_qa_string( &data[0], i, j ) != MB_SUCCESS ) return MB_FAILURE;
1767 
1768  qa_record_list.push_back( &data[0] );
1769  }
1770  }
1771 
1772  return MB_SUCCESS;
1773 }

References GET_DIM, max_str_length, MB_SUCCESS, and read_qa_string().

Referenced by read_qa_records().

◆ read_qa_records()

ErrorCode moab::ReadNCDF::read_qa_records ( EntityHandle  file_set)
private

Definition at line 1725 of file ReadNCDF.cpp.

1726 {
1727  std::vector< std::string > qa_records;
1728  read_qa_information( qa_records );
1729 
1730  std::vector< char > tag_data;
1731  for( std::vector< std::string >::iterator i = qa_records.begin(); i != qa_records.end(); ++i )
1732  {
1733  std::copy( i->begin(), i->end(), std::back_inserter( tag_data ) );
1734  tag_data.push_back( '\0' );
1735  }
1736 
1737  // If there were qa_records...tag them to the mCurrentMeshHandle
1738  if( !tag_data.empty() )
1739  {
1740  const void* ptr = &tag_data[0];
1741  int size = tag_data.size();
1742  if( mdbImpl->tag_set_by_ptr( mQaRecordTag, &file_set, 1, &ptr, &size ) != MB_SUCCESS )
1743  {
1744  return MB_FAILURE;
1745  }
1746  }
1747 
1748  return MB_SUCCESS;
1749 }

References MB_SUCCESS, mdbImpl, mQaRecordTag, read_qa_information(), size, and moab::Interface::tag_set_by_ptr().

Referenced by load_file().

◆ read_qa_string()

ErrorCode moab::ReadNCDF::read_qa_string ( char *  string,
int  record_number,
int  record_position 
)
private

Definition at line 1775 of file ReadNCDF.cpp.

1776 {
1777  std::vector< int > dims;
1778  int temp_var = -1;
1779  GET_VAR( "qa_records", temp_var, dims );
1780  if( -1 == temp_var )
1781  {
1782  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting qa record variable" );
1783  return MB_FAILURE;
1784  }
1785  size_t count[3], start[3];
1786  start[0] = record_number;
1787  start[1] = record_position;
1788  start[2] = 0;
1789  count[0] = 1;
1790  count[1] = 1;
1791  count[2] = max_str_length;
1792  int fail = nc_get_vara_text( ncFile, temp_var, start, count, temp_string );
1793  if( NC_NOERR != fail )
1794  {
1795  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem setting current record number variable position" );
1796  }
1797  // Get the variable id in the exodus file
1798 
1799  return MB_SUCCESS;
1800 }

References moab::fail(), GET_VAR, max_str_length, MB_SET_ERR, MB_SUCCESS, and ncFile.

Referenced by read_qa_information().

◆ read_sidesets()

ErrorCode moab::ReadNCDF::read_sidesets ( )
private

read the sidesets (does nothing for now)

Definition at line 1205 of file ReadNCDF.cpp.

1206 {
1207  // Uhhh if you read the same file (previously_read==true) then blow away the
1208  // sidesets pertaining to this file? How do you do that? If you read a
1209  // new file make sure all the offsets are correct.
1210 
1211  // If not loading any sidesets -- exit
1212  if( 0 == numberSideSets_loading ) return MB_SUCCESS;
1213 
1214  // Read in the sideset ids
1215  std::vector< int > id_array( numberSideSets_loading );
1216  int temp_var = -1;
1217  GET_1D_INT_VAR( "ss_prop1", temp_var, id_array );
1218  if( -1 == temp_var )
1219  {
1220  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting ss_prop1 variable" );
1221  }
1222 
1223  // Create a side set for each one
1224  int number_sides_in_set;
1225  int number_dist_factors_in_set;
1226 
1227  // Maybe there is already a sidesets meshset here we can append to
1228  Range child_meshsets;
1229  if( mdbImpl->get_entities_by_type( 0, MBENTITYSET, child_meshsets ) != MB_SUCCESS ) return MB_FAILURE;
1230 
1231  child_meshsets = subtract( child_meshsets, initRange );
1232 
1233  Range::iterator iter, end_iter;
1234 
1235  int i;
1236  std::vector< char > temp_string_storage( max_str_length + 1 );
1237  char* temp_string = &temp_string_storage[0];
1238  int temp_dim;
1239  for( i = 0; i < numberSideSets_loading; i++ )
1240  {
1241  // Get sideset parameters
1242  GET_DIMB( temp_dim, temp_string, "num_side_ss%d", i + 1, number_sides_in_set );
1243  GET_DIMB( temp_dim, temp_string, "num_df_ss%d", i + 1, number_dist_factors_in_set );
1244 
1245  // Size new arrays and get element and side lists
1246  std::vector< int > side_list( number_sides_in_set );
1247  std::vector< int > element_list( number_sides_in_set );
1248  INS_ID( temp_string, "side_ss%d", i + 1 );
1249  temp_var = -1;
1250  GET_1D_INT_VAR( temp_string, temp_var, side_list );
1251  if( -1 == temp_var )
1252  {
1253  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting sideset side variable" );
1254  }
1255 
1256  INS_ID( temp_string, "elem_ss%d", i + 1 );
1257  temp_var = -1;
1258  GET_1D_INT_VAR( temp_string, temp_var, element_list );
1259  if( -1 == temp_var )
1260  {
1261  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting sideset elem variable" );
1262  }
1263 
1264  std::vector< double > temp_dist_factor_vector;
1265  std::vector< EntityHandle > entities_to_add, reverse_entities;
1266  // Create the sideset entities
1267  if( create_ss_elements( &element_list[0], &side_list[0], number_sides_in_set, number_dist_factors_in_set,
1268  entities_to_add, reverse_entities, temp_dist_factor_vector, i + 1 ) != MB_SUCCESS )
1269  return MB_FAILURE;
1270 
1271  // If there are elements to add
1272  if( !entities_to_add.empty() || !reverse_entities.empty() )
1273  {
1274  iter = child_meshsets.begin();
1275  end_iter = child_meshsets.end();
1276 
1277  EntityHandle ss_handle = 0;
1278  for( ; iter != end_iter; ++iter )
1279  {
1280  int sideset_id;
1281  if( mdbImpl->tag_get_data( mNeumannSetTag, &( *iter ), 1, &sideset_id ) != MB_SUCCESS ) continue;
1282 
1283  if( id_array[i] == sideset_id )
1284  {
1285  // Found the meshset
1286  ss_handle = *iter;
1287  break;
1288  }
1289  }
1290 
1291  // If we didn't find a sideset already
1292  if( ss_handle == 0 )
1293  {
1294  if( mdbImpl->create_meshset( MESHSET_ORDERED | MESHSET_TRACK_OWNER, ss_handle ) != MB_SUCCESS )
1295  return MB_FAILURE;
1296 
1297  if( ss_handle == 0 ) return MB_FAILURE;
1298 
1299  int sideset_id = id_array[i];
1300  if( mdbImpl->tag_set_data( mNeumannSetTag, &ss_handle, 1, &sideset_id ) != MB_SUCCESS )
1301  return MB_FAILURE;
1302  if( mdbImpl->tag_set_data( mGlobalIdTag, &ss_handle, 1, &sideset_id ) != MB_SUCCESS ) return MB_FAILURE;
1303 
1304  if( !reverse_entities.empty() )
1305  {
1306  // Also make a reverse set to put in this set
1307  EntityHandle reverse_set;
1309  return MB_FAILURE;
1310 
1311  // Add the reverse set to the sideset set and the entities to the reverse set
1312  ErrorCode result = mdbImpl->add_entities( ss_handle, &reverse_set, 1 );
1313  if( MB_SUCCESS != result ) return result;
1314 
1315  result = mdbImpl->add_entities( reverse_set, &reverse_entities[0], reverse_entities.size() );
1316  if( MB_SUCCESS != result ) return result;
1317 
1318  // Set the reverse tag
1319  Tag sense_tag;
1320  int dum_sense = 0;
1321  result = mdbImpl->tag_get_handle( "NEUSET_SENSE", 1, MB_TYPE_INTEGER, sense_tag,
1322  MB_TAG_SPARSE | MB_TAG_CREAT, &dum_sense );
1323  if( result != MB_SUCCESS ) return result;
1324  dum_sense = -1;
1325  result = mdbImpl->tag_set_data( sense_tag, &reverse_set, 1, &dum_sense );
1326  if( result != MB_SUCCESS ) return result;
1327  }
1328  }
1329 
1330  if( mdbImpl->add_entities( ss_handle, ( entities_to_add.empty() ) ? NULL : &entities_to_add[0],
1331  entities_to_add.size() ) != MB_SUCCESS )
1332  return MB_FAILURE;
1333 
1334  // Distribution factor stuff
1335  if( number_dist_factors_in_set )
1336  {
1337  // If this sideset does not already has a distribution factor array...set one
1338  const void* ptr = 0;
1339  int size = 0;
1340  if( MB_SUCCESS == mdbImpl->tag_get_by_ptr( mDistFactorTag, &ss_handle, 1, &ptr, &size ) )
1341  {
1342  const double* data = reinterpret_cast< const double* >( ptr );
1343  std::copy( data, data + size, std::back_inserter( temp_dist_factor_vector ) );
1344  }
1345 
1346  ptr = &temp_dist_factor_vector[0];
1347  size = temp_dist_factor_vector.size();
1348  if( mdbImpl->tag_set_by_ptr( mDistFactorTag, &ss_handle, 1, &ptr, &size ) != MB_SUCCESS )
1349  return MB_FAILURE;
1350  }
1351  }
1352  }
1353 
1354  return MB_SUCCESS;
1355 }

References moab::Interface::add_entities(), moab::Range::begin(), moab::Interface::create_meshset(), create_ss_elements(), moab::Range::end(), ErrorCode, GET_1D_INT_VAR, GET_DIMB, moab::Interface::get_entities_by_type(), initRange, INS_ID, max_str_length, MB_SET_ERR, MB_SUCCESS, MB_TAG_CREAT, MB_TAG_SPARSE, MB_TYPE_INTEGER, MBENTITYSET, mdbImpl, mDistFactorTag, MESHSET_SET, MESHSET_TRACK_OWNER, mGlobalIdTag, mNeumannSetTag, numberSideSets_loading, size, moab::subtract(), moab::Interface::tag_get_by_ptr(), moab::Interface::tag_get_data(), moab::Interface::tag_get_handle(), moab::Interface::tag_set_by_ptr(), and moab::Interface::tag_set_data().

Referenced by load_file().

◆ read_tag_values()

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_nameThe file to read.
tag_nameThe tag for which to read values
tag_values_outOutput: The list of tag values.
subset_listAn 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_lengthThe length of the 'subset_list' array.

Implements moab::ReaderIface.

Definition at line 214 of file ReadNCDF.cpp.

219 {
220  if( subset_list )
221  {
222  MB_SET_ERR( MB_UNSUPPORTED_OPERATION, "ExodusII reader supports subset read only by material ID" );
223  }
224 
225  // Open netcdf/exodus file
226  int fail = nc_open( file_name, 0, &ncFile );
227  if( NC_NOWRITE != fail )
228  {
229  MB_SET_ERR( MB_FILE_DOES_NOT_EXIST, "ReadNCDF:: problem opening Netcdf/Exodus II file " << file_name );
230  }
231 
232  // 1. Read the header
233  ErrorCode rval = read_exodus_header();
234  if( MB_FAILURE == rval ) return rval;
235 
236  int count = 0;
237  const char* prop;
238  const char* blocks = "eb_prop1";
239  const char* nodesets = "ns_prop1";
240  const char* sidesets = "ss_prop1";
241 
242  if( !strcmp( tag_name, MATERIAL_SET_TAG_NAME ) )
243  {
245  prop = blocks;
246  }
247  else if( !strcmp( tag_name, DIRICHLET_SET_TAG_NAME ) )
248  {
249  count = numberNodeSets_loading;
250  prop = nodesets;
251  }
252  else if( !strcmp( tag_name, NEUMANN_SET_TAG_NAME ) )
253  {
254  count = numberSideSets_loading;
255  prop = sidesets;
256  }
257  else
258  {
259  ncFile = 0;
260  return MB_TAG_NOT_FOUND;
261  }
262 
263  if( count )
264  {
265  int nc_var = -1;
266  GET_1D_INT_VAR( prop, nc_var, id_array );
267  if( !nc_var )
268  {
269  MB_SET_ERR( MB_FAILURE, "Problem getting prop variable" );
270  }
271  }
272 
273  // Close the file
274  fail = nc_close( ncFile );
275  if( NC_NOERR != fail )
276  {
277  MB_SET_ERR( MB_FAILURE, "Trouble closing file" );
278  }
279 
280  ncFile = 0;
281  return MB_SUCCESS;
282 }

References DIRICHLET_SET_TAG_NAME, ErrorCode, moab::fail(), GET_1D_INT_VAR, MATERIAL_SET_TAG_NAME, MB_FILE_DOES_NOT_EXIST, MB_SET_ERR, MB_SUCCESS, MB_TAG_NOT_FOUND, MB_UNSUPPORTED_OPERATION, ncFile, NEUMANN_SET_TAG_NAME, numberElementBlocks_loading, numberNodeSets_loading, numberSideSets_loading, and read_exodus_header().

◆ reset()

◆ tokenize()

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

Definition at line 2448 of file ReadNCDF.cpp.

2449 {
2450  std::string::size_type last = str.find_first_not_of( delimiters, 0 );
2451  std::string::size_type pos = str.find_first_of( delimiters, last );
2452  while( std::string::npos != pos && std::string::npos != last )
2453  {
2454  tokens.push_back( str.substr( last, pos - last ) );
2455  last = str.find_first_not_of( delimiters, pos );
2456  pos = str.find_first_of( delimiters, last );
2457  if( std::string::npos == pos ) pos = str.size();
2458  }
2459 }

Referenced by update().

◆ update()

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 1805 of file ReadNCDF.cpp.

1810 {
1811  // Function : updating current database from new exodus_file.
1812  // Creator: Jane Hu
1813  // opts is currently designed as following
1814  // tdata = <var_name>[, time][,op][,destination]
1815  // where var_name show the tag name to be updated, this version just takes
1816  // coord.
1817  // time is the optional, and it gives time step of each of the mesh
1818  // info in exodus file. It start from 1.
1819  // op is the operation that is going to be performed on the var_name info.
1820  // currently support 'set'
1821  // destination shows where to store the updated info, currently assume it is
1822  // stored in the same database by replacing the old info if there's no input
1823  // for destination, or the destination data is given in exodus format and just
1824  // need to update the coordinates.
1825  // Assumptions:
1826  // 1. Assume the num_el_blk's in both DB and update exodus file are the same.
1827  // 2. Assume num_el_in_blk1...num_el_in_blk(num_el_blk) numbers are matching, may in
1828  // different order. example: num_el_in_blk11 = num_el_in_blk22 && num_el_in_blk12 =
1829  // num_el_in_blk21.
1830  // 3. In exodus file, get node_num_map
1831  // 4. loop through the node_num_map, use it to find the node in the cub file.
1832  // 5. Replace coord[0][n] with coordx[m] + vals_nod_var1(time_step, m) for all directions for
1833  // matching nodes.
1834  // Test: try to get hexes
1835 
1836  // *******************************************************************
1837  // Move nodes to their deformed locations.
1838  // *******************************************************************
1839  ErrorCode rval;
1840  std::string s;
1841  rval = opts.get_str_option( "tdata", s );
1842  if( MB_SUCCESS != rval )
1843  {
1844  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem reading file options" );
1845  }
1846  std::vector< std::string > tokens;
1847  tokenize( s, tokens, "," );
1848 
1849  // 1. Check for time step to find the match time
1850  int time_step = 1;
1851  if( tokens.size() > 1 && !tokens[1].empty() )
1852  {
1853  const char* time_s = tokens[1].c_str();
1854  char* endptr;
1855  long int pval = strtol( time_s, &endptr, 0 );
1856  std::string end = endptr;
1857  if( !end.empty() ) // Syntax error
1858  return MB_TYPE_OUT_OF_RANGE;
1859 
1860  // Check for overflow (parsing long int, returning int)
1861  time_step = pval;
1862  if( pval != (long int)time_step ) return MB_TYPE_OUT_OF_RANGE;
1863  if( time_step <= 0 ) return MB_TYPE_OUT_OF_RANGE;
1864  }
1865 
1866  // 2. Check for the operations, currently support set.
1867  const char* op;
1868  if( tokens.size() < 3 || ( tokens[2] != "set" && tokens[2] != "add" ) )
1869  {
1870  MB_SET_ERR( MB_TYPE_OUT_OF_RANGE, "ReadNCDF: invalid operation specified for update" );
1871  }
1872  op = tokens[2].c_str();
1873 
1874  // Open netcdf/exodus file
1875  ncFile = 0;
1876  int fail = nc_open( exodus_file_name, 0, &ncFile );
1877  if( !ncFile )
1878  {
1879  MB_SET_ERR( MB_FILE_DOES_NOT_EXIST, "ReadNCDF:: problem opening Netcdf/Exodus II file " << exodus_file_name );
1880  }
1881 
1882  rval = read_exodus_header();MB_CHK_ERR( rval );
1883 
1884  // Check to make sure that the requested time step exists
1885  int ncdim = -1;
1886  int max_time_steps;
1887  GET_DIM( ncdim, "time_step", max_time_steps );
1888  if( -1 == ncdim )
1889  {
1890  std::cout << "ReadNCDF: could not get number of time steps" << std::endl;
1891  return MB_FAILURE;
1892  }
1893  std::cout << " Maximum time step=" << max_time_steps << std::endl;
1894  if( max_time_steps < time_step )
1895  {
1896  std::cout << "ReadNCDF: time step is greater than max_time_steps" << std::endl;
1897  return MB_FAILURE;
1898  }
1899 
1900  // Get the time
1901  std::vector< double > times( max_time_steps );
1902  int nc_var = -1;
1903  GET_1D_DBL_VAR( "time_whole", nc_var, times );
1904  if( -1 == nc_var )
1905  {
1906  std::cout << "ReadNCDF: unable to get time variable" << std::endl;
1907  }
1908  else
1909  {
1910  std::cout << " Step " << time_step << " is at " << times[time_step - 1] << " seconds" << std::endl;
1911  }
1912 
1913  // Read in the node_num_map.
1914  std::vector< int > ptr( numberNodes_loading );
1915 
1916  int varid = -1;
1917  GET_1D_INT_VAR( "node_num_map", varid, ptr );
1918  if( -1 == varid )
1919  {
1920  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting node number map data" );
1921  }
1922 
1923  // Read in the deformations.
1924  std::vector< std::vector< double > > deformed_arrays( 3 );
1925  std::vector< std::vector< double > > orig_coords( 3 );
1926  deformed_arrays[0].reserve( numberNodes_loading );
1927  deformed_arrays[1].reserve( numberNodes_loading );
1928  deformed_arrays[2].reserve( numberNodes_loading );
1929  orig_coords[0].reserve( numberNodes_loading );
1930  orig_coords[1].reserve( numberNodes_loading );
1931  orig_coords[2].reserve( numberNodes_loading );
1932  size_t start[2] = { static_cast< size_t >( time_step - 1 ), 0 },
1933  count[2] = { 1, static_cast< size_t >( numberNodes_loading ) };
1934  std::vector< int > dims;
1935  int coordx = -1, coordy = -1, coordz = -1;
1936  GET_VAR( "vals_nod_var1", coordx, dims );
1937  GET_VAR( "vals_nod_var2", coordy, dims );
1938  if( numberDimensions_loading == 3 ) GET_VAR( "vals_nod_var3", coordz, dims );
1939  if( -1 == coordx || -1 == coordy || ( numberDimensions_loading == 3 && -1 == coordz ) )
1940  {
1941  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting coords variable" );
1942  }
1943 
1944  fail = nc_get_vara_double( ncFile, coordx, start, count, &deformed_arrays[0][0] );
1945  if( NC_NOERR != fail )
1946  {
1947  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting x deformation array" );
1948  }
1949 
1950  fail = nc_get_vara_double( ncFile, coordy, start, count, &deformed_arrays[1][0] );
1951  if( NC_NOERR != fail )
1952  {
1953  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting y deformation array" );
1954  }
1955  if( numberDimensions_loading == 3 )
1956  {
1957  fail = nc_get_vara_double( ncFile, coordz, start, count, &deformed_arrays[2][0] );
1958  if( NC_NOERR != fail )
1959  {
1960  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting z deformation array" );
1961  }
1962  }
1963 
1964  int coord1 = -1, coord2 = -1, coord3 = -1;
1965  GET_1D_DBL_VAR( "coordx", coord1, orig_coords[0] );
1966  if( -1 == coord1 )
1967  {
1968  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting x coord array" );
1969  }
1970  GET_1D_DBL_VAR( "coordy", coord2, orig_coords[1] );
1971  if( -1 == coord2 )
1972  {
1973  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting y coord array" );
1974  }
1975  if( numberDimensions_loading == 3 )
1976  {
1977  GET_1D_DBL_VAR( "coordz", coord3, orig_coords[2] );
1978  if( -1 == coord3 )
1979  {
1980  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting z coord array" );
1981  }
1982  }
1983 
1984  // b. Deal with DB file : get node info. according to node_num_map.
1985  if( tokens[0] != "coord" && tokens[0] != "COORD" ) return MB_NOT_IMPLEMENTED;
1986 
1987  if( strcmp( op, "set" ) != 0 && strcmp( op, " set" ) != 0 ) return MB_NOT_IMPLEMENTED;
1988 
1989  // Two methods of matching nodes (id vs. proximity)
1990  const bool match_node_ids = true;
1991 
1992  // Get nodes in cubit file
1993  Range cub_verts;
1994  rval = mdbImpl->get_entities_by_type( cub_file_set, MBVERTEX, cub_verts );MB_CHK_ERR( rval );
1995  std::cout << " cub_file_set contains " << cub_verts.size() << " nodes." << std::endl;
1996 
1997  // Some accounting
1998  std::cout << " exodus file contains " << numberNodes_loading << " nodes." << std::endl;
1999  double max_magnitude = 0;
2000  double average_magnitude = 0;
2001  int found = 0;
2002  int lost = 0;
2003  std::map< int, EntityHandle > cub_verts_id_map;
2004  AdaptiveKDTree kdtree( mdbImpl );
2005  EntityHandle root;
2006 
2007  // Should not use cub verts unless they have been matched. Place in a map
2008  // for fast handle_by_id lookup.
2009  std::map< int, EntityHandle > matched_cub_vert_id_map;
2010 
2011  // Place cub verts in a map for searching by id
2012  if( match_node_ids )
2013  {
2014  std::vector< int > cub_ids( cub_verts.size() );
2015  rval = mdbImpl->tag_get_data( mGlobalIdTag, cub_verts, &cub_ids[0] );MB_CHK_ERR( rval );
2016  for( unsigned i = 0; i != cub_verts.size(); ++i )
2017  {
2018  cub_verts_id_map.insert( std::pair< int, EntityHandle >( cub_ids[i], cub_verts[i] ) );
2019  }
2020 
2021  // Place cub verts in a kdtree for searching by proximity
2022  }
2023  else
2024  {
2025  FileOptions tree_opts( "MAX_PER_LEAF=1;SPLITS_PER_DIR=1;CANDIDATE_PLANE_SET=0" );
2026  rval = kdtree.build_tree( cub_verts, &root, &tree_opts );MB_CHK_ERR( rval );
2027  AdaptiveKDTreeIter tree_iter;
2028  rval = kdtree.get_tree_iterator( root, tree_iter );MB_CHK_ERR( rval );
2029  }
2030 
2031  // For each exo vert, find the matching cub vert
2032  for( int i = 0; i < numberNodes_loading; ++i )
2033  {
2034  int exo_id = ptr[i];
2035  CartVect exo_coords( orig_coords[0][i], orig_coords[1][i], orig_coords[2][i] );
2036  EntityHandle cub_vert = 0;
2037  bool found_match = false;
2038 
2039  // By id
2040  if( match_node_ids )
2041  {
2042  std::map< int, EntityHandle >::iterator i_iter;
2043  i_iter = cub_verts_id_map.find( exo_id );
2044  if( i_iter != cub_verts_id_map.end() )
2045  {
2046  found_match = true;
2047  cub_vert = i_iter->second;
2048  }
2049 
2050  // By proximity
2051  }
2052  else
2053  {
2054  // The MAX_NODE_DIST is the farthest distance to search for a node.
2055  // For the 1/12th symmetry 85 pin model, the max node dist could not be less
2056  // than 1e-1 (March 26, 2010).
2057  const double MAX_NODE_DIST = 1e-1;
2058 
2059  std::vector< EntityHandle > leaves;
2060  double min_dist = MAX_NODE_DIST;
2061  rval = kdtree.distance_search( exo_coords.array(), MAX_NODE_DIST, leaves );MB_CHK_ERR( rval );
2062  for( std::vector< EntityHandle >::const_iterator j = leaves.begin(); j != leaves.end(); ++j )
2063  {
2064  std::vector< EntityHandle > leaf_verts;
2065  rval = mdbImpl->get_entities_by_type( *j, MBVERTEX, leaf_verts );MB_CHK_ERR( rval );
2066  for( std::vector< EntityHandle >::const_iterator k = leaf_verts.begin(); k != leaf_verts.end(); ++k )
2067  {
2068  CartVect orig_cub_coords, difference;
2069  rval = mdbImpl->get_coords( &( *k ), 1, orig_cub_coords.array() );MB_CHK_ERR( rval );
2070  difference = orig_cub_coords - exo_coords;
2071  double dist = difference.length();
2072  if( dist < min_dist )
2073  {
2074  min_dist = dist;
2075  cub_vert = *k;
2076  }
2077  }
2078  }
2079  if( 0 != cub_vert ) found_match = true;
2080  }
2081 
2082  // If a match is found, update it with the deformed coords from the exo file.
2083  if( found_match )
2084  {
2085  CartVect updated_exo_coords;
2086  matched_cub_vert_id_map.insert( std::pair< int, EntityHandle >( exo_id, cub_vert ) );
2087  updated_exo_coords[0] = orig_coords[0][i] + deformed_arrays[0][i];
2088  updated_exo_coords[1] = orig_coords[1][i] + deformed_arrays[1][i];
2089 
2090  if( numberDimensions_loading == 3 ) updated_exo_coords[2] = orig_coords[2][i] + deformed_arrays[2][i];
2091  rval = mdbImpl->set_coords( &cub_vert, 1, updated_exo_coords.array() );MB_CHK_ERR( rval );
2092  ++found;
2093 
2094  double magnitude =
2095  sqrt( deformed_arrays[0][i] * deformed_arrays[0][i] + deformed_arrays[1][i] * deformed_arrays[1][i] +
2096  deformed_arrays[2][i] * deformed_arrays[2][i] );
2097  if( magnitude > max_magnitude ) max_magnitude = magnitude;
2098  average_magnitude += magnitude;
2099  }
2100  else
2101  {
2102  ++lost;
2103  std::cout << "cannot match exo node " << exo_id << " " << exo_coords << std::endl;
2104  }
2105  }
2106 
2107  // Exo verts that could not be matched have already been printed. Now print the
2108  // cub verts that could not be matched.
2109  if( matched_cub_vert_id_map.size() < cub_verts.size() )
2110  {
2111  Range unmatched_cub_verts = cub_verts;
2112  for( std::map< int, EntityHandle >::const_iterator i = matched_cub_vert_id_map.begin();
2113  i != matched_cub_vert_id_map.end(); ++i )
2114  {
2115  unmatched_cub_verts.erase( i->second );
2116  }
2117 
2118  for( Range::const_iterator i = unmatched_cub_verts.begin(); i != unmatched_cub_verts.end(); ++i )
2119  {
2120  int cub_id;
2121  rval = mdbImpl->tag_get_data( mGlobalIdTag, &( *i ), 1, &cub_id );MB_CHK_ERR( rval );
2122 
2123  CartVect cub_coords;
2124  rval = mdbImpl->get_coords( &( *i ), 1, cub_coords.array() );MB_CHK_ERR( rval );
2125  std::cout << "cannot match cub node " << cub_id << " " << cub_coords << std::endl;
2126  }
2127 
2128  std::cout << " " << unmatched_cub_verts.size() << " nodes from the cub file could not be matched."
2129  << std::endl;
2130  }
2131 
2132  // Summarize statistics
2133  std::cout << " " << found << " nodes from the exodus file were matched in the cub_file_set ";
2134  if( match_node_ids )
2135  {
2136  std::cout << "by id." << std::endl;
2137  }
2138  else
2139  {
2140  std::cout << "by proximity." << std::endl;
2141  }
2142 
2143  // Fail if all of the nodes could not be matched.
2144  if( 0 != lost )
2145  {
2146  std::cout << "Error: " << lost << " nodes from the exodus file could not be matched." << std::endl;
2147  // return MB_FAILURE;
2148  }
2149  std::cout << " maximum node displacement magnitude: " << max_magnitude << " cm" << std::endl;
2150  std::cout << " average node displacement magnitude: " << average_magnitude / found << " cm" << std::endl;
2151 
2152  // *******************************************************************
2153  // Remove dead elements from the MOAB instance.
2154  // *******************************************************************
2155 
2156  // How many element variables are in the file?
2157  int n_elem_var;
2158  GET_DIM( ncdim, "num_elem_var", n_elem_var );
2159 
2160  // Get element variable names
2161  varid = -1;
2162  int cstatus = nc_inq_varid( ncFile, "name_elem_var", &varid );
2163  std::vector< char > names_memory( n_elem_var * max_str_length );
2164  std::vector< char* > names( n_elem_var );
2165  for( int i = 0; i < n_elem_var; ++i )
2166  names[i] = &names_memory[i * max_str_length];
2167  if( cstatus != NC_NOERR || varid == -1 )
2168  {
2169  std::cout << "ReadNCDF: name_elem_var does not exist" << std::endl;
2170  return MB_FAILURE;
2171  }
2172  int status = nc_get_var_text( ncFile, varid, &names_memory[0] );
2173  if( NC_NOERR != status )
2174  {
2175  MB_SET_ERR( MB_FAILURE, "ReadNCDF: Problem getting element variable names" );
2176  }
2177 
2178  // Is one of the element variable names "death_status"? If so, get its index
2179  // in the element variable array.
2180  int death_index;
2181  bool found_death_index = false;
2182  for( int i = 0; i < n_elem_var; ++i )
2183  {
2184  std::string temp( names[i] );
2185  std::string::size_type pos0 = temp.find( "death_status" );
2186  std::string::size_type pos1 = temp.find( "Death_Status" );
2187  std::string::size_type pos2 = temp.find( "DEATH_STATUS" );
2188  if( std::string::npos != pos0 || std::string::npos != pos1 || std::string::npos != pos2 )
2189  {
2190  found_death_index = true;
2191  death_index = i + 1; // NetCDF variables start with 1
2192  break;
2193  }
2194  }
2195  if( !found_death_index )
2196  {
2197  MB_SET_ERR( MB_FAILURE, "ReadNCDF: Problem getting index of death_status variable" );
2198  }
2199 
2200  // The exodus header has already been read. This contains the number of element
2201  // blocks.
2202 
2203  // Dead elements are listed by block. Read the block headers to determine how
2204  // many elements are in each block.
2205  rval = read_block_headers( blocks_to_load, num_blocks );
2206  if( MB_FAILURE == rval )
2207  {
2208  MB_SET_ERR( MB_FAILURE, "ReadNCDF: Problem reading block headers" );
2209  }
2210 
2211  // Dead elements from the Exodus file can be located in the cub_file_set by id
2212  // or by connectivity. Currently, finding elements by id requires careful book
2213  // keeping when constructing the model in Cubit. To avoid this, one can match
2214  // elements by connectivity instead.
2215  const bool match_elems_by_connectivity = true;
2216 
2217  // Get the element id map. The ids in the map are from the elements in the blocks.
2218  // elem_num_map(blk1 elem ids, blk2 elem ids, blk3 elem ids, ...)
2219  std::vector< int > elem_ids( numberNodes_loading );
2220  if( !match_elems_by_connectivity )
2221  {
2222  GET_1D_INT_VAR( "elem_num_map", varid, elem_ids );
2223  if( -1 == varid )
2224  {
2225  MB_SET_ERR( MB_FAILURE, "ReadNCDF: Problem getting element number map data" );
2226  }
2227  }
2228 
2229  // For each block
2230  int first_elem_id_in_block = 0;
2231  int block_count = 1; // NetCDF variables start with 1
2232  int total_elems = 0;
2233  int total_dead_elems = 0;
2234  for( std::vector< ReadBlockData >::iterator i = blocksLoading.begin(); i != blocksLoading.end(); ++i )
2235  {
2236  // Get the ncdf connect variable
2237  std::string temp_string( "connect" );
2238  std::stringstream temp_ss;
2239  temp_ss << block_count;
2240  temp_string += temp_ss.str();
2241  temp_string += "\0";
2242  std::vector< int > ldims;
2243  GET_VAR( temp_string.c_str(), nc_var, ldims );
2244  if( !nc_var )
2245  {
2246  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting connect variable" );
2247  }
2248  // The element type is an attribute of the connectivity variable
2249  nc_type att_type;
2250  size_t att_len;
2251  fail = nc_inq_att( ncFile, nc_var, "elem_type", &att_type, &att_len );
2252  if( NC_NOERR != fail )
2253  {
2254  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting elem type attribute" );
2255  }
2256  std::vector< char > dum_str( att_len + 1 );
2257  fail = nc_get_att_text( ncFile, nc_var, "elem_type", &dum_str[0] );
2258  if( NC_NOERR != fail )
2259  {
2260  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting elem type" );
2261  }
2262  ExoIIElementType elem_type = ExoIIUtil::static_element_name_to_type( &dum_str[0] );
2263  ( *i ).elemType = elem_type;
2264  const EntityType mb_type = ExoIIUtil::ExoIIElementMBEntity[( *i ).elemType];
2265 
2266  // Get the number of nodes per element
2267  unsigned int nodes_per_element = ExoIIUtil::VerticesPerElement[( *i ).elemType];
2268 
2269  // Read the connectivity into that memory.
2270  // int exo_conn[i->numElements][nodes_per_element];
2271  int* exo_conn = new int[i->numElements * nodes_per_element];
2272  // NcBool status = temp_var->get(&exo_conn[0][0], i->numElements, nodes_per_element);
2273  size_t lstart[2] = { 0, 0 }, lcount[2];
2274  lcount[0] = i->numElements;
2275  lcount[1] = nodes_per_element;
2276  fail = nc_get_vara_int( ncFile, nc_var, lstart, lcount, exo_conn );
2277  if( NC_NOERR != fail )
2278  {
2279  delete[] exo_conn;
2280  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting connectivity" );
2281  }
2282 
2283  // Get the death_status at the correct time step.
2284  std::vector< double > death_status( i->numElements ); // It seems wrong, but it uses doubles
2285  std::string array_name( "vals_elem_var" );
2286  temp_ss.str( "" ); // stringstream won't clear by temp.clear()
2287  temp_ss << death_index;
2288  array_name += temp_ss.str();
2289  array_name += "eb";
2290  temp_ss.str( "" ); // stringstream won't clear by temp.clear()
2291  temp_ss << block_count;
2292  array_name += temp_ss.str();
2293  array_name += "\0";
2294  GET_VAR( array_name.c_str(), nc_var, ldims );
2295  if( !nc_var )
2296  {
2297  MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting death status variable" );
2298  }
2299  lstart[0] = time_step - 1;
2300  lstart[1] = 0;
2301  lcount[0] = 1;
2302  lcount[1] = i->numElements;
2303  status = nc_get_vara_double( ncFile, nc_var, lstart, lcount, &death_status[0] );
2304  if( NC_NOERR != status )
2305  {
2306  delete[] exo_conn;
2307  MB_SET_ERR( MB_FAILURE, "ReadNCDF: Problem setting time step for death_status" );
2308  }
2309 
2310  // Look for dead elements. If there is too many dead elements and this starts
2311  // to take too long, I should put the elems in a kd-tree for more efficient
2312  // searching. Alternatively I could get the exo connectivity and match nodes.
2313  int dead_elem_counter = 0, missing_elem_counter = 0;
2314  for( int j = 0; j < i->numElements; ++j )
2315  {
2316  if( 1 != death_status[j] )
2317  {
2318  Range cub_elem, cub_nodes;
2319  if( match_elems_by_connectivity )
2320  {
2321  // Get exodus nodes for the element
2322  std::vector< int > elem_conn( nodes_per_element );
2323  for( unsigned int k = 0; k < nodes_per_element; ++k )
2324  {
2325  // elem_conn[k] = exo_conn[j][k];
2326  elem_conn[k] = exo_conn[j * nodes_per_element + k];
2327  }
2328  // Get the ids of the nodes (assume we are matching by id)
2329  // Remember that the exodus array locations start with 1 (not 0).
2330  std::vector< int > elem_conn_node_ids( nodes_per_element );
2331  for( unsigned int k = 0; k < nodes_per_element; ++k )
2332  {
2333  elem_conn_node_ids[k] = ptr[elem_conn[k] - 1];
2334  }
2335  // Get the cub_file_set nodes by id
2336  // The map is a log search and takes almost no time.
2337  // MOAB's linear tag search takes 5-10 minutes.
2338  for( unsigned int k = 0; k < nodes_per_element; ++k )
2339  {
2340  std::map< int, EntityHandle >::iterator k_iter;
2341  k_iter = matched_cub_vert_id_map.find( elem_conn_node_ids[k] );
2342 
2343  if( k_iter == matched_cub_vert_id_map.end() )
2344  {
2345  std::cout << "ReadNCDF: Found no cub node with id=" << elem_conn_node_ids[k]
2346  << ", but expected to find only 1." << std::endl;
2347  break;
2348  }
2349  cub_nodes.insert( k_iter->second );
2350  }
2351 
2352  if( nodes_per_element != cub_nodes.size() )
2353  {
2354  std::cout << "ReadNCDF: nodes_per_elemenet != cub_nodes.size()" << std::endl;
2355  delete[] exo_conn;
2356  return MB_INVALID_SIZE;
2357  }
2358 
2359  // Get the cub_file_set element with the same nodes
2360  int to_dim = CN::Dimension( mb_type );
2361  rval = mdbImpl->get_adjacencies( cub_nodes, to_dim, false, cub_elem );
2362  if( MB_SUCCESS != rval )
2363  {
2364  std::cout << "ReadNCDF: could not get dead cub element" << std::endl;
2365  delete[] exo_conn;
2366  return rval;
2367  }
2368 
2369  // Pronto/Presto renumbers elements, so matching cub and exo elements by
2370  // id is not possible at this time.
2371  }
2372  else
2373  {
2374  // Get dead element's id
2375  int elem_id = elem_ids[first_elem_id_in_block + j];
2376  void* id[] = { &elem_id };
2377  // Get the element by id
2378  rval = mdbImpl->get_entities_by_type_and_tag( cub_file_set, mb_type, &mGlobalIdTag, id, 1, cub_elem,
2380  if( MB_SUCCESS != rval )
2381  {
2382  delete[] exo_conn;
2383  return rval;
2384  }
2385  }
2386 
2387  if( 1 == cub_elem.size() )
2388  {
2389  // Delete the dead element from the cub file. It will be removed from sets
2390  // ONLY if they are tracking meshsets.
2391  rval = mdbImpl->remove_entities( cub_file_set, cub_elem );
2392  if( MB_SUCCESS != rval )
2393  {
2394  delete[] exo_conn;
2395  return rval;
2396  }
2397  rval = mdbImpl->delete_entities( cub_elem );
2398  if( MB_SUCCESS != rval )
2399  {
2400  delete[] exo_conn;
2401  return rval;
2402  }
2403  }
2404  else
2405  {
2406  std::cout << "ReadNCDF: Should have found 1 element with type=" << mb_type
2407  << " in cub_file_set, but instead found " << cub_elem.size() << std::endl;
2408  rval = mdbImpl->list_entities( cub_nodes );
2409  ++missing_elem_counter;
2410  delete[] exo_conn;
2411  return MB_FAILURE;
2412  }
2413  ++dead_elem_counter;
2414  }
2415  }
2416  // Print some statistics
2417  temp_ss.str( "" );
2418  temp_ss << i->blockId;
2419  total_dead_elems += dead_elem_counter;
2420  total_elems += i->numElements;
2421  std::cout << " Block " << temp_ss.str() << " has " << dead_elem_counter << "/" << i->numElements
2422  << " dead elements." << std::endl;
2423  if( 0 != missing_elem_counter )
2424  {
2425  std::cout << " " << missing_elem_counter
2426  << " dead elements in this block were not found in the cub_file_set." << std::endl;
2427  }
2428 
2429  // Advance the pointers into element ids and block_count. Memory cleanup.
2430  first_elem_id_in_block += i->numElements;
2431  ++block_count;
2432  delete[] exo_conn;
2433  }
2434 
2435  std::cout << " Total: " << total_dead_elems << "/" << total_elems << " dead elements." << std::endl;
2436 
2437  // Close the file
2438  fail = nc_close( ncFile );
2439  if( NC_NOERR != fail )
2440  {
2441  MB_SET_ERR( MB_FAILURE, "Trouble closing file" );
2442  }
2443 
2444  ncFile = 0;
2445  return MB_SUCCESS;
2446 }

References moab::CartVect::array(), moab::Range::begin(), blocksLoading, moab::AdaptiveKDTree::build_tree(), moab::Interface::delete_entities(), moab::CN::Dimension(), moab::AdaptiveKDTree::distance_search(), moab::Range::end(), moab::Range::erase(), ErrorCode, moab::ExoIIUtil::ExoIIElementMBEntity, moab::fail(), GET_1D_DBL_VAR, GET_1D_INT_VAR, moab::Interface::get_adjacencies(), moab::Interface::get_coords(), GET_DIM, moab::Interface::get_entities_by_type(), moab::Interface::get_entities_by_type_and_tag(), moab::FileOptions::get_str_option(), moab::AdaptiveKDTree::get_tree_iterator(), GET_VAR, moab::Range::insert(), moab::Interface::INTERSECT, moab::CartVect::length(), moab::Interface::list_entities(), max_str_length, MB_CHK_ERR, MB_FILE_DOES_NOT_EXIST, MB_INVALID_SIZE, MB_NOT_IMPLEMENTED, MB_SET_ERR, MB_SUCCESS, MB_TYPE_OUT_OF_RANGE, MBVERTEX, mdbImpl, mGlobalIdTag, ncFile, numberDimensions_loading, numberNodes_loading, read_block_headers(), read_exodus_header(), moab::Interface::remove_entities(), moab::Interface::set_coords(), moab::Range::size(), moab::ExoIIUtil::static_element_name_to_type(), moab::Interface::tag_get_data(), tokenize(), and moab::ExoIIUtil::VerticesPerElement.

Referenced by load_file().

Member Data Documentation

◆ blocksLoading

std::vector< ReadBlockData > moab::ReadNCDF::blocksLoading
private

◆ CPU_WORD_SIZE

int moab::ReadNCDF::CPU_WORD_SIZE
private

Definition at line 194 of file ReadNCDF.hpp.

Referenced by read_exodus_header().

◆ exodusFile

std::string moab::ReadNCDF::exodusFile
private

file name

Definition at line 201 of file ReadNCDF.hpp.

◆ faceBlocksLoading

std::vector< ReadFaceBlockData > moab::ReadNCDF::faceBlocksLoading
private

Definition at line 236 of file ReadNCDF.hpp.

Referenced by read_face_blocks_headers(), read_polyhedra_faces(), and reset().

◆ initRange

Range moab::ReadNCDF::initRange
private

range of entities in initial mesh, before this read

Definition at line 251 of file ReadNCDF.hpp.

Referenced by load_file(), read_nodesets(), and read_sidesets().

◆ IO_WORD_SIZE

int moab::ReadNCDF::IO_WORD_SIZE
private

Definition at line 195 of file ReadNCDF.hpp.

Referenced by read_exodus_header().

◆ max_line_length

int moab::ReadNCDF::max_line_length
private

Definition at line 248 of file ReadNCDF.hpp.

Referenced by read_exodus_header().

◆ max_str_length

◆ mCurrentMeshHandle

EntityHandle moab::ReadNCDF::mCurrentMeshHandle
private

Meshset Handle for the mesh that is currently being read.

Definition at line 224 of file ReadNCDF.hpp.

Referenced by reset().

◆ mdbImpl

◆ mDirichletSetTag

Tag moab::ReadNCDF::mDirichletSetTag
private

Definition at line 241 of file ReadNCDF.hpp.

Referenced by read_nodesets(), and ReadNCDF().

◆ mDistFactorTag

Tag moab::ReadNCDF::mDistFactorTag
private

Definition at line 244 of file ReadNCDF.hpp.

Referenced by read_nodesets(), read_sidesets(), and ReadNCDF().

◆ mGlobalIdTag

Tag moab::ReadNCDF::mGlobalIdTag
private

◆ mHasMidNodesTag

Tag moab::ReadNCDF::mHasMidNodesTag
private

Definition at line 243 of file ReadNCDF.hpp.

Referenced by read_elements(), and ReadNCDF().

◆ mMaterialSetTag

Tag moab::ReadNCDF::mMaterialSetTag
private

Cached tags for reading. Note that all these tags are defined when the core is initialized.

Definition at line 240 of file ReadNCDF.hpp.

Referenced by read_elements(), and ReadNCDF().

◆ mNeumannSetTag

Tag moab::ReadNCDF::mNeumannSetTag
private

Definition at line 242 of file ReadNCDF.hpp.

Referenced by read_sidesets(), and ReadNCDF().

◆ mQaRecordTag

Tag moab::ReadNCDF::mQaRecordTag
private

Definition at line 246 of file ReadNCDF.hpp.

Referenced by read_qa_records(), and ReadNCDF().

◆ ncFile

int moab::ReadNCDF::ncFile
private

◆ nodesInLoadedBlocks

std::vector< char > moab::ReadNCDF::nodesInLoadedBlocks
private

Definition at line 227 of file ReadNCDF.hpp.

Referenced by read_elements(), read_nodesets(), read_polyhedra_faces(), and reset().

◆ numberDimensions_loading

int moab::ReadNCDF::numberDimensions_loading
private

Definition at line 221 of file ReadNCDF.hpp.

Referenced by number_dimensions(), read_exodus_header(), read_nodes(), reset(), and update().

◆ numberElementBlocks_loading

int moab::ReadNCDF::numberElementBlocks_loading
private

number of blocks in the current exoII file

Definition at line 210 of file ReadNCDF.hpp.

Referenced by read_block_headers(), read_exodus_header(), read_tag_values(), and reset().

◆ numberElements_loading

int moab::ReadNCDF::numberElements_loading
private

number of elements in the current exoII file

Definition at line 207 of file ReadNCDF.hpp.

Referenced by read_exodus_header(), read_global_ids(), and reset().

◆ numberFaceBlocks_loading

int moab::ReadNCDF::numberFaceBlocks_loading
private

number of face blocks in the current exoII file (used for polyhedra)

Definition at line 213 of file ReadNCDF.hpp.

Referenced by load_file(), read_exodus_header(), read_face_blocks_headers(), and reset().

◆ numberNodes_loading

int moab::ReadNCDF::numberNodes_loading
private

number of nodes in the current exoII file

Definition at line 204 of file ReadNCDF.hpp.

Referenced by read_elements(), read_exodus_header(), read_global_ids(), read_nodes(), read_polyhedra_faces(), reset(), and update().

◆ numberNodeSets_loading

int moab::ReadNCDF::numberNodeSets_loading
private

number of nodesets in the current exoII file

Definition at line 216 of file ReadNCDF.hpp.

Referenced by read_exodus_header(), read_nodesets(), read_tag_values(), and reset().

◆ numberSideSets_loading

int moab::ReadNCDF::numberSideSets_loading
private

number of sidesets in the current exoII file

Definition at line 219 of file ReadNCDF.hpp.

Referenced by read_exodus_header(), read_sidesets(), read_tag_values(), and reset().

◆ polyfaces

std::vector< EntityHandle > moab::ReadNCDF::polyfaces
private

Definition at line 233 of file ReadNCDF.hpp.

Referenced by read_elements(), and read_polyhedra_faces().

◆ readMeshIface

ReadUtilIface* moab::ReadNCDF::readMeshIface
private

Definition at line 108 of file ReadNCDF.hpp.

Referenced by read_elements(), read_nodes(), ReadNCDF(), and ~ReadNCDF().

◆ vertexOffset

int moab::ReadNCDF::vertexOffset
private

int to offset vertex ids with

Definition at line 198 of file ReadNCDF.hpp.

Referenced by read_elements(), read_global_ids(), read_nodes(), read_nodesets(), read_polyhedra_faces(), and reset().


The documentation for this class was generated from the following files: