Commit db5298c9 authored by Kenneth Leiter's avatar Kenneth Leiter
Browse files

ENH: Add XdmfExodusReader for reading ExodusII files into Xdmf. Add ability...

ENH: Add XdmfExodusReader for reading ExodusII files into Xdmf.  Add ability for XdmfArray::initialize() to take a size value to initialize the array to.
parent 732cc4b3
#
# Find the Exodus finite element data model library from Sandia
#
# EXODUS_FOUND - System has Exodus
# EXODUS_INCLUDE_DIR - The LibXml2 include directory
# EXODUS_LIBRARIES - The libraries needed to use LibXml2
FIND_PACKAGE(NetCDF REQUIRED)
FIND_PATH(EXODUS_INCLUDE_DIR NAMES exodusII.h)
FIND_LIBRARY(EXODUS_LIBRARIES NAMES exodusii exodusIIv2c)
INCLUDE(FindPackageHandleStandardArgs)
# handle the QUIETLY and REQUIRED arguments and set EXODUS_FOUND to TRUE if
# all listed variables are TRUE
FIND_PACKAGE_HANDLE_STANDARD_ARGS(Exodus DEFAULT_MSG EXODUS_LIBRARIES EXODUS_INCLUDE_DIR)
MARK_AS_ADVANCED(EXODUS_INCLUDE_DIR EXODUS_LIBRARIES)
#
# Find NetCDF include directories and libraries
#
# NetCDF_FOUND - System has NetCDF
# NetCDF_INCLUDE_DIR - The NetCDF include directory
# NetCDF_LIBRARIES - The libraries needed to use NetCDF
FIND_PATH(NetCDF_INCLUDE_DIR netcdf.h)
FIND_LIBRARY(NetCDF_LIBRARIES
NAMES netcdf
${NetCDF_PREFIX}
${NetCDF_PREFIX}/lib64
${NetCDF_PREFIX}/lib
/usr/local/lib64
/usr/lib64
/usr/lib64/netcdf-3
/usr/local/lib
/usr/lib
/usr/lib/netcdf-3
)
INCLUDE(FindPackageHandleStandardArgs)
# handle the QUIETLY and REQUIRED arguments and set EXODUS_FOUND to TRUE if
# all listed variables are TRUE
FIND_PACKAGE_HANDLE_STANDARD_ARGS(NetCDF DEFAULT_MSG NetCDF_INCLUDE_DIR NetCDF_LIBRARIES)
MARK_AS_ADVANCED(NetCDF_INCLUDE_DIR NetCDF_LIBRARIES)
......@@ -95,6 +95,12 @@ boost::shared_ptr<const XdmfTopologyType> XdmfTopologyType::Quadrilateral_8()
return p;
}
boost::shared_ptr<const XdmfTopologyType> XdmfTopologyType::Quadrilateral_9()
{
static boost::shared_ptr<const XdmfTopologyType> p(new XdmfTopologyType(9, "Quadrilateral_9", Quadratic));
return p;
}
boost::shared_ptr<const XdmfTopologyType> XdmfTopologyType::Tetrahedron_10()
{
static boost::shared_ptr<const XdmfTopologyType> p(new XdmfTopologyType(10, "Tetrahedron_10", Quadratic));
......@@ -113,6 +119,12 @@ boost::shared_ptr<const XdmfTopologyType> XdmfTopologyType::Wedge_15()
return p;
}
boost::shared_ptr<const XdmfTopologyType> XdmfTopologyType::Wedge_18()
{
static boost::shared_ptr<const XdmfTopologyType> p(new XdmfTopologyType(18, "Wedge_18", Quadratic));
return p;
}
boost::shared_ptr<const XdmfTopologyType> XdmfTopologyType::Hexahedron_20()
{
static boost::shared_ptr<const XdmfTopologyType> p(new XdmfTopologyType(20, "Hexahedron_20", Quadratic));
......
......@@ -23,6 +23,7 @@
* Edge_3
* Triangle_6
* Quadrilateral_8
* Quadrilateral_9
* Tetrahedron_10
* Pyramid_13
* Wedge_15
......@@ -37,7 +38,7 @@
* TwoDCoRectMesh
* ThreeDSMesh
* ThreeDRectMesh
* ThreeDCoRectMesh
* ThreeDCoRectMesh
*/
class XdmfTopologyType : public XdmfItemProperty {
......@@ -65,9 +66,11 @@ public:
static boost::shared_ptr<const XdmfTopologyType> Edge_3();
static boost::shared_ptr<const XdmfTopologyType> Triangle_6();
static boost::shared_ptr<const XdmfTopologyType> Quadrilateral_8();
static boost::shared_ptr<const XdmfTopologyType> Quadrilateral_9();
static boost::shared_ptr<const XdmfTopologyType> Tetrahedron_10();
static boost::shared_ptr<const XdmfTopologyType> Pyramid_13();
static boost::shared_ptr<const XdmfTopologyType> Wedge_15();
static boost::shared_ptr<const XdmfTopologyType> Wedge_18();
static boost::shared_ptr<const XdmfTopologyType> Hexahedron_20();
static boost::shared_ptr<const XdmfTopologyType> Hexahedron_24();
static boost::shared_ptr<const XdmfTopologyType> Hexahedron_27();
......
......@@ -471,43 +471,43 @@ std::string XdmfArray::getValuesString() const
return "";
}
void XdmfArray::initialize(const boost::shared_ptr<const XdmfArrayType> arrayType)
void XdmfArray::initialize(const boost::shared_ptr<const XdmfArrayType> arrayType, const unsigned int size)
{
if(arrayType == XdmfArrayType::Int8())
{
this->initialize<char>();
this->initialize<char>(size);
}
else if(arrayType == XdmfArrayType::Int16())
{
this->initialize<short>();
this->initialize<short>(size);
}
else if(arrayType == XdmfArrayType::Int32())
{
this->initialize<int>();
this->initialize<int>(size);
}
else if(arrayType == XdmfArrayType::Int64())
{
this->initialize<long>();
this->initialize<long>(size);
}
else if(arrayType == XdmfArrayType::Float32())
{
this->initialize<float>();
this->initialize<float>(size);
}
else if(arrayType == XdmfArrayType::Float64())
{
this->initialize<double>();
this->initialize<double>(size);
}
else if(arrayType == XdmfArrayType::UInt8())
{
this->initialize<unsigned char>();
this->initialize<unsigned char>(size);
}
else if(arrayType == XdmfArrayType::UInt16())
{
this->initialize<unsigned short>();
this->initialize<unsigned short>(size);
}
else if(arrayType == XdmfArrayType::UInt32())
{
this->initialize<unsigned int>();
this->initialize<unsigned int>(size);
}
else if(arrayType == XdmfArrayType::Uninitialized())
{
......
......@@ -187,15 +187,19 @@ public:
/**
* Initializes the array to contain an empty container of a particular type.
*
* @param size the number of values in the initialized array.
* @return a smart pointer to the internal vector of values initialized in this array.
*/
template <typename T>
boost::shared_ptr<std::vector<T> > initialize();
boost::shared_ptr<std::vector<T> > initialize(const unsigned int size = 0);
/**
* Initializes the array to contain an empty container of a particular XdmfArrayType.
*
* @param arrayType the type of array to initialize.
* @param size the number of values in the initialized array.
*/
void initialize(const boost::shared_ptr<const XdmfArrayType> arrayType);
void initialize(const boost::shared_ptr<const XdmfArrayType> arrayType, const unsigned int size = 0);
/**
* Returns whether the array is initialized (contains values in memory).
......
......@@ -198,14 +198,14 @@ void XdmfArray::getValuesCopy(const unsigned int startIndex, T * const valuesPoi
}
template <typename T>
boost::shared_ptr<std::vector<T> > XdmfArray::initialize()
boost::shared_ptr<std::vector<T> > XdmfArray::initialize(const unsigned int size)
{
if(mHaveArrayPointer)
{
releaseArrayPointer();
}
// Set type of variant to type of pointer
boost::shared_ptr<std::vector<T> > newArray(new std::vector<T>());
boost::shared_ptr<std::vector<T> > newArray(new std::vector<T>(size));
if(mTmpReserveSize > 0)
{
newArray->reserve(mTmpReserveSize);
......
......@@ -49,8 +49,7 @@ void XdmfHDF5Controller::read(XdmfArray * const array)
hssize_t numVals = H5Sget_simple_extent_npoints(dataspace);
hid_t datatype = H5Dget_type(dataset);
array->initialize(mType);
array->resize(numVals, 0);
array->initialize(mType, numVals);
H5Dread(dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, array->getValuesPointer());
herr_t status;
......
include_directories(${CMAKE_CURRENT_SOURCE_DIR})
option(XDMF_BUILD_EXODUS_IO OFF)
option(XDMF_BUILD_PARTITIONER OFF)
set(XdmfUtilsExecutables)
......@@ -8,6 +9,15 @@ set(XdmfUtilsSources
)
set(XdmfUtilsLinkLibraries Xdmf)
if(XDMF_BUILD_EXODUS_IO)
find_package(Exodus REQUIRED)
if(EXODUS_FOUND)
include_directories(${EXODUS_INCLUDE_DIR} ${NetCDF_INCLUDE_DIR})
endif(EXODUS_FOUND)
set(XdmfUtilsSources ${XdmfUtilsSources} XdmfExodusReader)
set(XdmfUtilsLinkLibraries ${XdmfUtilsLinkLibraries} ${EXODUS_LIBRARIES} ${NetCDF_LIBRARIES})
endif(XDMF_BUILD_EXODUS_IO)
if(XDMF_BUILD_PARTITIONER)
find_package(Metis REQUIRED)
if(METIS_FOUND)
......
/*******************************************************************/
/* XDMF */
/* eXtensible Data Model and Format */
/* */
/* Id : Id */
/* Date : $Date$ */
/* Version : $Revision$ */
/* */
/* Author: */
/* Kenneth Leiter */
/* kenneth.leiter@arl.army.mil */
/* US Army Research Laboratory */
/* Aberdeen Proving Ground, MD */
/* */
/* Copyright @ 2009 US Army Research Laboratory */
/* All Rights Reserved */
/* See Copyright.txt or http://www.arl.hpc.mil/ice for details */
/* */
/* This software is distributed WITHOUT ANY WARRANTY; without */
/* even the implied warranty of MERCHANTABILITY or FITNESS */
/* FOR A PARTICULAR PURPOSE. See the above copyright notice */
/* for more information. */
/* */
/*******************************************************************/
#include <exodusII.h>
#include "XdmfArrayType.hpp"
#include "XdmfAttribute.hpp"
#include "XdmfAttributeCenter.hpp"
#include "XdmfAttributeType.hpp"
#include "XdmfExodusReader.hpp"
#include "XdmfGeometry.hpp"
#include "XdmfGeometryType.hpp"
#include "XdmfGrid.hpp"
#include "XdmfSet.hpp"
#include "XdmfSetType.hpp"
#include "XdmfTopology.hpp"
#include "XdmfTopologyType.hpp"
boost::shared_ptr<XdmfExodusReader> XdmfExodusReader::New()
{
boost::shared_ptr<XdmfExodusReader> p(new XdmfExodusReader());
return p;
}
XdmfExodusReader::XdmfExodusReader()
{
}
XdmfExodusReader::~XdmfExodusReader()
{
}
boost::shared_ptr<XdmfGrid> XdmfExodusReader::read(const std::string & fileName) const
{
boost::shared_ptr<XdmfGrid> grid = XdmfGrid::New();
// Read Exodus II file to XdmfGrid via Exodus II API
float version;
int CPU_word_size = sizeof(double);
int IO_word_size = 0; // Get from file
int exodusHandle = ex_open(fileName.c_str(), EX_READ, &CPU_word_size, &IO_word_size, &version);
char * title = new char[MAX_LINE_LENGTH+1];
int num_dim, num_nodes, num_elem, num_elem_blk, num_node_sets, num_side_sets;
ex_get_init (exodusHandle, title, &num_dim, &num_nodes, &num_elem, &num_elem_blk, &num_node_sets, &num_side_sets);
delete [] title;
/*
cout << "Title: " << title <<
"\nNum Dim: " << num_dim <<
"\nNum Nodes: " << num_nodes <<
"\nNum Elem: " << num_elem <<
"\nNum Elem Blk: " << num_elem_blk <<
"\nNum Node Sets: " << num_node_sets <<
"\nNum Side Sets: " << num_side_sets << endl;
*/
// Read geometry values
double * x = new double[num_nodes];
double * y = new double[num_nodes];
double * z = new double[num_nodes];
ex_get_coord(exodusHandle, x, y, z);
// In the future we may want to do XDMF_GEOMETRY_X_Y_Z?
if(num_dim == 3)
{
grid->getGeometry()->setType(XdmfGeometryType::XYZ());
}
else if(num_dim = 2)
{
grid->getGeometry()->setType(XdmfGeometryType::XY());
}
else
{
// Xdmf does not support geometries with less than 2 dimensions
assert(false);
}
grid->getGeometry()->initialize(XdmfArrayType::Float64());
grid->getGeometry()->reserve(num_nodes * num_dim);
for(unsigned int i=0; i<num_nodes; ++i)
{
grid->getGeometry()->pushBack(x[i]);
grid->getGeometry()->pushBack(y[i]);
if(num_dim == 3)
{
grid->getGeometry()->pushBack(z[i]);
}
}
delete [] x;
delete [] y;
delete [] z;
int * blockIds = new int[num_elem_blk];
ex_get_elem_blk_ids(exodusHandle, blockIds);
int * numElemsInBlock = new int[num_elem_blk];
int * numNodesPerElemInBlock = new int[num_elem_blk];
int * numElemAttrInBlock = new int[num_elem_blk];
std::vector<boost::shared_ptr<const XdmfTopologyType> > topologyTypes;
topologyTypes.reserve(num_elem_blk);
int totalNumElem = 0;
int totalConns = 0;
for(unsigned int i=0; i<num_elem_blk; ++i)
{
char * elem_type = new char[MAX_STR_LENGTH+1];
int num_nodes_per_elem, num_elem_this_blk, num_attr;
ex_get_elem_block(exodusHandle, blockIds[i], elem_type, &num_elem_this_blk, &num_nodes_per_elem, &num_attr);
/*
cout << "Block Id: " << blockIds[j] <<
"\nElem Type: " << elem_type <<
"\nNum Elem in Blk: " << num_elem_this_blk <<
"\nNum Nodes per Elem: " << num_nodes_per_elem <<
"\nNum Attr: " << num_attr << endl;
*/
numElemsInBlock[i] = num_elem_this_blk;
numNodesPerElemInBlock[i] = num_nodes_per_elem;
numElemAttrInBlock[i] = num_attr;
topologyTypes.push_back(this->exodusToXdmfTopologyType(elem_type, num_nodes_per_elem));
totalNumElem += num_elem_this_blk;
totalConns += num_elem_this_blk * num_nodes_per_elem;
delete [] elem_type;
}
if(topologyTypes.size() > 0)
{
grid->getTopology()->setType(topologyTypes[0]);
if(topologyTypes.size() > 1)
{
for(std::vector<boost::shared_ptr<const XdmfTopologyType> >::const_iterator iter = topologyTypes.begin() + 1; iter != topologyTypes.end(); ++iter)
{
// Cannot be mixed topology!
assert(grid->getTopology()->getType() == *iter);
}
}
}
topologyTypes.clear();
grid->getTopology()->initialize(XdmfArrayType::Int32(), totalConns);
int * connectivityPointer = (int *)grid->getTopology()->getValuesPointer();
// Read connectivity from element blocks
int elemIndex = 0;
for(unsigned int i=0; i<num_elem_blk; ++i)
{
ex_get_elem_conn(exodusHandle, blockIds[i], &connectivityPointer[elemIndex]);
elemIndex += numElemsInBlock[i] * numNodesPerElemInBlock[i];
}
// This is taken from VTK's vtkExodusIIReader and adapted to fit Xdmf element types, which have the same ordering as VTK.
if(grid->getTopology()->getType() == XdmfTopologyType::Hexahedron_20() || grid->getTopology()->getType() == XdmfTopologyType::Hexahedron_27())
{
int * ptr = connectivityPointer;
int itmp[4];
// Exodus Node ordering does not match Xdmf, we must convert.
for(unsigned int i=0; i<totalNumElem; ++i)
{
ptr += 12;
for(unsigned int j=0; j<4; ++j, ++ptr)
{
itmp[j] = *ptr;
*ptr = ptr[4];
}
for(unsigned int j=0; j<4; ++j, ++ptr)
{
*ptr = itmp[j];
}
if(grid->getTopology()->getType() == XdmfTopologyType::Hexahedron_27())
{
for(unsigned int j=0; j<4; ++j, ++ptr)
{
itmp[j] = *ptr;
*ptr = ptr[3];
}
*(ptr++) = itmp[1];
*(ptr++) = itmp[2];
*(ptr++) = itmp[0];
}
}
}
else if(grid->getTopology()->getType() == XdmfTopologyType::Wedge_15() || grid->getTopology()->getType() == XdmfTopologyType::Wedge_18())
{
int * ptr = connectivityPointer;
int itmp[3];
// Exodus Node ordering does not match Xdmf, we must convert.
for (unsigned int i=0; i<totalNumElem; i++)
{
ptr += 9;
for(unsigned int j=0; j<3; ++j, ++ptr)
{
itmp[j] = *ptr;
*ptr = ptr[3];
}
for(unsigned int j=0; j<3; ++j, ++ptr)
{
*ptr = itmp[j];
}
if(grid->getTopology()->getType() == XdmfTopologyType::Wedge_18())
{
itmp[0] = *(ptr);
itmp[1] = *(ptr+1);
itmp[2] = *(ptr+2);
*(ptr++) = itmp[1];
*(ptr++) = itmp[2];
*(ptr++) = itmp[0];
}
}
}
// Subtract all node ids by 1 since exodus indices start at 1
for(unsigned int i=0; i<totalConns; ++i)
{
connectivityPointer[i]--;
}
boost::shared_ptr<XdmfAttribute> globalIds = XdmfAttribute::New();
globalIds->setName("GlobalNodeId");
globalIds->setCenter(XdmfAttributeCenter::Node());
globalIds->setType(XdmfAttributeType::GlobalId());
globalIds->initialize(XdmfArrayType::Int32(), num_nodes);
int * globalIdsPointer = (int*)globalIds->getValuesPointer();
ex_get_node_num_map(exodusHandle, globalIdsPointer);
// Subtract all node ids by 1 since exodus indices start at 1
for(unsigned int i=0; i<num_nodes; ++i)
{
globalIdsPointer[i]--;
}
grid->insert(globalIds);
// Read node sets
int * nodeSetIds = new int[num_node_sets];
ex_get_node_set_ids(exodusHandle, nodeSetIds);
char * node_set_names[num_node_sets];
for (unsigned int i=0; i<num_node_sets; ++i)
{
node_set_names[i] = new char[MAX_STR_LENGTH+1];
}
ex_get_names(exodusHandle, EX_NODE_SET, node_set_names);
for (unsigned int i=0; i<num_node_sets; ++i)
{
int num_nodes_in_set, num_df_in_set;
ex_get_node_set_param(exodusHandle, nodeSetIds[i], &num_nodes_in_set, &num_df_in_set);
/*
cout << "Node Set Id: " << nodeSetIds[j] <<
"\nNode Set Name: " << node_set_names[j] <<
"\nNum Nodes in Set: "<< num_nodes_in_set <<
"\nNum Distrub Factors: " << num_df_in_set << endl;
*/
if (num_nodes_in_set > 0)
{
int * node_set_node_list = new int[num_nodes_in_set];
ex_get_node_set(exodusHandle, nodeSetIds[i], node_set_node_list);
boost::shared_ptr<XdmfSet> set = XdmfSet::New();
set->setName(node_set_names[i]);
set->setType(XdmfSetType::Node());
for(unsigned int j=0; j<num_nodes_in_set; ++j)
{
set->insert(node_set_node_list[j] - 1);
}
grid->insert(set);
/*
if(num_df_in_set > 0)
{
double * node_set_distribution_factors = new double[num_df_in_set];
ex_get_node_set_dist_fact(exodusHandle, nodeSetIds[j], node_set_distribution_factors);
XdmfAttribute * attr = new XdmfAttribute();
attr->SetName("SetAttribute");
attr->SetAttributeType(XDMF_ATTRIBUTE_TYPE_SCALAR);
attr->SetAttributeCenter(XDMF_ATTRIBUTE_CENTER_NODE);
attr->SetDeleteOnGridDelete(true);
XdmfArray * attrVals = attr->GetValues();
attrVals->SetNumberType(XDMF_FLOAT32_TYPE);
attrVals->SetNumberOfElements(num_df_in_set);
attrVals->SetValues(0, node_set_distribution_factors, num_df_in_set, 1, 1);
set->Insert(attr);
delete [] node_set_distribution_factors;
}
*/
delete [] node_set_node_list;
}
delete [] node_set_names[i];
}
delete [] nodeSetIds;
// Read result variables (attributes)
int num_global_vars, num_nodal_vars, num_elem_vars;
ex_get_var_param(exodusHandle, "g", &num_global_vars);
ex_get_var_param(exodusHandle, "n", &num_nodal_vars);
ex_get_var_param(exodusHandle, "e", &num_elem_vars);
/*
cout << "Num Global Vars: " << num_global_vars <<
"\nNum Nodal Vars: " << num_nodal_vars <<
"\nNum Elem Vars: " << num_elem_vars << endl;
*/
char * global_var_names[num_global_vars];
char * nodal_var_names[num_nodal_vars];
char * elem_var_names[num_elem_vars];
for (int j=0; j<num_global_vars; j++)
{
global_var_names[j] = new char[MAX_STR_LENGTH+1];
}
for (int j=0; j<num_nodal_vars; j++)
{
nodal_var_names[j] = new char[MAX_STR_LENGTH+1];
}
for (int j=0; j<num_elem_vars; j++)
{
elem_var_names[j] = new char[MAX_STR_LENGTH+1];
}
ex_get_var_names(exodusHandle, "g", num_global_vars, global_var_names);
ex_get_var_names(exodusHandle, "n", num_nodal_vars, nodal_var_names);
ex_get_var_names(exodusHandle, "e", num_elem_vars, elem_var_names);
/*
cout << "Global Vars Names: " << endl;
for (int j=0; j<num_global_vars; j++)
{
cout << global_var_names[j] << endl;
}
cout << "Nodal Vars Names: " << endl;
for (int j=0; j<num_nodal_vars; j++)
{
cout << nodal_var_names[j] << endl;
}
cout << "Elem Vars Names: " << endl;