Mesh Oriented datABase  (version 5.5.1)
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, max_str_length + 1 );
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, max_str_length + 1 );
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, max_str_length + 1 );
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, max_str_length + 1 );
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, max_str_length + 1 );
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, max_str_length + 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, max_str_length + 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, max_str_length + 1 );
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, max_str_length + 1 );
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, max_str_length + 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, max_str_length + 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 2445 of file ReadNCDF.cpp.

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

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

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

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: