[cig-commits] r13838 - cs/cigma/trunk/src
luis at geodynamics.org
luis at geodynamics.org
Mon Jan 12 14:10:06 PST 2009
Author: luis
Date: 2009-01-12 14:10:06 -0800 (Mon, 12 Jan 2009)
New Revision: 13838
Added:
cs/cigma/trunk/src/io_netcdf.cpp
cs/cigma/trunk/src/io_netcdf.h
Modified:
cs/cigma/trunk/src/io_exo_reader.cpp
cs/cigma/trunk/src/io_exo_reader.h
Log:
Added ExodusReader class for handling CUBIT meshes
For now, instead of using the ExodusII API, we directly use the
NetCDF API (isolating these functions in the io_netcdf.{h,cpp}
files).
Note that we create a cigma::MeshPart object by using the first
element block in the file (this needs to be generalized to using
multiple blocks, possibly each having a different element type).
Modified: cs/cigma/trunk/src/io_exo_reader.cpp
===================================================================
--- cs/cigma/trunk/src/io_exo_reader.cpp 2009-01-12 22:10:04 UTC (rev 13837)
+++ cs/cigma/trunk/src/io_exo_reader.cpp 2009-01-12 22:10:06 UTC (rev 13838)
@@ -1,125 +1,292 @@
-#include "ExodusII_Reader.h"
+#include "io_exo_reader.h"
+#include "io_netcdf.h"
+#ifdef HAVE_NETCDF
+//#include <algorithm>
+//#include <cctype>
+#include <cstring>
+
+#include "nc_array.h"
+#include "eb_array.h"
+
+
+using namespace std;
using namespace cigma;
+using boost::shared_ptr;
-ExodusII_Reader::ExodusII_Reader()
+ExodusReader::ExodusReader()
{
}
-ExodusII_Reader::~ExodusII_Reader()
+ExodusReader::~ExodusReader()
{
}
-int ExodusII_Reader::open(const char *filename, const char *mode)
+int ExodusReader::open(const char *filename, const char *mode)
{
- TRI_LOG_STR("ExodusII_Reader::open()");
+ TRI_LOG_STR("ExodusReader::open()");
TRI_LOG(filename);
- TRI_LOG(mode);
+ //TRI_LOG(mode);
+ // XXX: use mode to determined second argument to NcFile constructor
this->filename = filename;
+ this->ncfile = shared_ptr<NcFile>(new NcFile(filename, NcFile::ReadOnly));
- FooBar reader
- if (reader)
+ if (!(ncfile->is_valid()))
{
- this->filename = filename;
- this->reader->Update();
- this->dataset = reader->GetOutput();
+ return -1;
}
+
+ return 0;
}
-int ExodusII_Reader::close()
+int ExodusReader::close()
{
}
-int ExodusII_Reader::getDataset(const char *loc, double **data, int *num, int *dim)
+int ExodusReader::getDataset(const char *loc, double **data, int *num, int *dim)
{
- TRI_LOG_STR("ExodusII_Reader::getDataset()");
+ TRI_LOG_STR("ExodusReader::getDataset()");
+ throw cigma::Exception("ExodusReader::getDataset", "Not implemented");
}
-int ExodusII_Reader::getCoordinates(const char *loc, double **coordinates, int *nno, int *nsd)
+int ExodusReader::getCoordinates(const char *loc, double **coordinates, int *nno, int *nsd)
{
- TRI_LOG_STR("ExodusII_Reader::getCoordinates()");
+ TRI_LOG_STR("ExodusReader::getCoordinates()");
+ valarray<double> coords;
+ nc_read_coordinates(*ncfile, coords, *nno, *nsd);
+ int size = (*nno) * (*nsd);
+ *coordinates = new double[size];
+ std::memcpy(*coordinates, &coords[0], size * (sizeof(double)));
+ /* XXX: debug output (remove)
+ cout << "nno = " << (*nno) << endl;
+ cout << "nsd = " << (*nsd) << endl;
+ for (int i = 0; i < size; i++)
+ {
+ (*coordinates)[i] = coords[i];
+ } // */
+ /*
+ for (int i = 0; i < (*nno); i++)
+ {
+ for (int j = 0; j < (*nsd); j++)
+ {
+ int n = (*nsd)*i + j;
+ cout << " " << coords[n];
+ }
+ cout << endl;
+ } // */
+ return 0;
}
-int ExodusII_Reader::getConnectivity(const char *loc, int **connectivity, int *nel, int *ndofs)
+int ExodusReader::getConnectivity(const char *loc, int **connectivity, int *nel, int *ndofs)
{
- TRI_LOG_STR("ExodusII_Reader::getConnectivity()");
+ TRI_LOG_STR("ExodusReader::getConnectivity()");
+ valarray<int> cells;
+ nc_read_connectivity(*ncfile, cells, *nel, *ndofs);
+ int size = (*nel) * (*ndofs);
+ *connectivity = new int[size];
+ std::memcpy(*connectivity, &cells[0], size * sizeof(int));
+ /* XXX: debug output (remove)
+ for (int i = 0; i < size; i++)
+ {
+ (*connectivity)[i] = cells[i];
+ } // */
+ /*
+ cout << "nel = " << (*nel) << endl;
+ cout << "ndofs = " << (*ndofs) << endl;
+ for (int i = 0; i < (*nel); i++)
+ {
+ for (int j = 0; j < (*ndofs); j++)
+ {
+ int n = (*ndofs)*i + j;
+ cout << " " << cells[n];
+ }
+ cout << endl;
+ } // */
+ return 0;
}
+void ExodusReader::getIntArray(const char *loc, cigma::array<int>& A)
+{
+ TRI_LOG_STR("ExodusReader::getIntArray()");
+ throw cigma::Exception("ExodusReader::getIntArray", "Not implemented");
+}
-void ExodusII_Reader::getIntArray(const char *loc, cigma::array<int>& A)
+
+void ExodusReader::getLongArray(const char *loc, cigma::array<long>& A)
{
- TRI_LOG_STR("ExodusII_Reader::getIntArray()");
+ TRI_LOG_STR("ExodusReader::getLongArray()");
+ throw cigma::Exception("ExodusReader::getLongArray", "Not implemented");
}
-void ExodusII_Reader::getLongArray(const char *loc, cigma::array<long>& A)
+void ExodusReader::getFloatArray(const char *loc, cigma::array<float>& A)
{
- TRI_LOG_STR("ExodusII_Reader::getIntArray()");
+ TRI_LOG_STR("ExodusReader::getFloatArray()");
+ throw cigma::Exception("ExodusReader::getFloatArray", "Not implemented");
}
-void ExodusII_Reader::getFloatArray(const char *loc, cigma::array<float>& A)
+void ExodusReader::getDoubleArray(const char *loc, cigma::array<double>& A)
{
- TRI_LOG_STR("ExodusII_Reader::getIntArray()");
+ TRI_LOG_STR("ExodusReader::getDoubleArray()");
+ throw cigma::Exception("ExodusReader::getDoubleArray", "Not implemented");
}
-void ExodusII_Reader::getDoubleArray(const char *loc, cigma::array<double>& A)
+/*
+static int lower_case(int c)
{
- TRI_LOG_STR("ExodusII_Reader::getIntArray()");
+ return tolower(c);
}
+*/
+Cell::type ExodusReader::getCellType(const char *loc)
+{
+ TRI_LOG_STR("ExodusReader::getCellType()");
-shared_ptr<NodeCoordinates> ExodusII_Reader::getNodeCoordinates(const char *loc)
+ string cellname, celltype;
+ nc_read_elem_type(*ncfile, loc, cellname);
+ TRI_LOG(cellname);
+
+ // use lowercase form of cellname
+ //transform(cellname.begin(), cellname.end(), cellname.begin(), lower_case);
+ //return Cell::string2type(cellname);
+
+ // alternatively, consider each case separately
+ if (cellname == "HEX8")
+ {
+ return Cell::HEX8;
+ }
+ else if (cellname == "TETRA")
+ {
+ return Cell::TET4;
+ }
+ return Cell::NONE;
+}
+
+
+shared_ptr<NodeCoordinates> ExodusReader::getNodeCoordinates(const char *loc)
{
- shared_ptr<nc_array> nc;
+ TRI_LOG_STR("ExodusReader::getNodeCoordinates()");
+ shared_ptr<nc_array> nc(new nc_array);
+ try
+ {
+ int status = this->getCoordinates(loc, &(nc->coords), &(nc->nno), &(nc->nsd));
+ if (status < 0)
+ {
+ std::ostringstream stream;
+ stream << "Could not read coordinates from exodus file '"
+ << filename << "'" << std::ends;
+ throw cigma::Exception("ExodusReader::getNodeCoordinates", stream.str());
+ }
+ }
+ catch (runtime_error& e)
+ {
+ std::ostringstream stream;
+ stream << "In exodus file '" << filename << "': "
+ << e.what() << std::ends;
+ throw cigma::Exception("ExodusReader::getNodeCoordinates", stream.str());
+ }
return nc;
}
-shared_ptr<ElementBlock> ExodusII_Reader::getElementBlock(const char *loc)
+shared_ptr<ElementBlock> ExodusReader::getElementBlock(const char *loc)
{
- shared_ptr<eb_array> eb;
+ TRI_LOG_STR("ExodusReader::getElementBlock()");
+ shared_ptr<eb_array> eb(new eb_array);
+ try
+ {
+ int status = this->getConnectivity(loc, &(eb->connect), &(eb->nel), &(eb->ndofs));
+ if (status < 0)
+ {
+ std::ostringstream stream;
+ stream << "Could not read coordinates from exodus file '"
+ << filename << "'" << std::ends;
+ throw cigma::Exception("ExodusReader::getNodeCoordinates", stream.str());
+ }
+ }
+ catch (runtime_error& e)
+ {
+ std::ostringstream stream;
+ stream << "In exodus file '" << filename << "': "
+ << e.what() << std::ends;
+ throw cigma::Exception("ExodusReader::getElementBlock", stream.str());
+ }
return eb;
}
-shared_ptr<MeshPart> ExodusII_Reader::getMeshPart(const char *loc)
+shared_ptr<MeshPart> ExodusReader::getMeshPart(const char *loc)
{
- shared_ptr<MeshPart> mesh;
+ shared_ptr<MeshPart> mesh(new MeshPart);
+ mesh->coords = this->getNodeCoordinates("coords");
+ mesh->connect = this->getElementBlock("connect1");
+ mesh->cell_type = this->getCellType("connect1");
return mesh;
}
-shared_ptr<Quadrature> ExodusII_Reader::getQuadrature(const char *loc)
+shared_ptr<Quadrature> ExodusReader::getQuadrature(const char *loc)
{
+ throw cigma::Exception("ExodusReader::getQuadrature", "Not implemented");
shared_ptr<Quadrature> Q;
return Q;
}
-shared_ptr<FE> ExodusII_Reader::getFE(const char *loc)
+shared_ptr<FE> ExodusReader::getFE(const char *loc)
{
+ throw cigma::Exception("ExodusReader::getFE", "Not implemented");
shared_ptr<FE> fe;
return fe;
}
-shared_ptr<DofHandler> ExodusII_Reader::getDofHandler(const char *loc)
+shared_ptr<DofHandler> ExodusReader::getDofHandler(const char *loc)
{
+ throw cigma::Exception("ExodusReader::getDofHandler", "Not implemented");
shared_ptr<DofHandler> dh;
return dh;
}
+int exo_print_contents(const char *filename)
+{
+ throw cigma::Exception("exo_print_contents", "Not implemented");
+ return 0;
+}
+
+#else
+
+
+cigma::ExodusReader::ExodusReader()
+{
+ std::string msg("Can't create ExodusReader. Cigma is not configured with the NetCDF C++ library");
+ throw cigma::Exception("ExodusReader", msg);
+}
+
+cigma::ExodusReader::~ExodusReader()
+{
+}
+
+int exo_print_contents(const char *filename)
+{
+ throw cigma::Exception("exo_print_contents",
+ "Error: Cannot list Exodus file. "
+ "Cigma not configured with NetCDF library", 1);
+ return -1;
+}
+
+#endif
Modified: cs/cigma/trunk/src/io_exo_reader.h
===================================================================
--- cs/cigma/trunk/src/io_exo_reader.h 2009-01-12 22:10:04 UTC (rev 13837)
+++ cs/cigma/trunk/src/io_exo_reader.h 2009-01-12 22:10:06 UTC (rev 13838)
@@ -1,24 +1,25 @@
-#ifndef __CIGMA_EXODUSII_READER_H__
-#define __CIGMA_EXODUSII_READER_H__
+#ifndef __CIGMA_EXODUS_READER_H__
+#define __CIGMA_EXODUS_READER_H__
#include "io_file_reader.h"
namespace cigma
{
- class ExodusII_Reader;
+ class ExodusReader;
}
+// XXX: move this function elsewhere?
+int exo_print_contents(const char *filename);
-#ifdef HAVE_VTK
+#ifdef HAVE_NETCDF
-#include "vtkExodusIIReader.h"
+#include <netcdfcpp.h>
-
-class cigma::ExodusII_Reader : public cigma::FileReader
+class cigma::ExodusReader : public cigma::FileReader
{
public:
- ExodusII_Reader();
- ~ExodusII_Reader();
+ ExodusReader();
+ ~ExodusReader();
ReaderType getReaderType() { return EXO_FILE_READER; }
@@ -34,6 +35,8 @@
void getFloatArray(const char *loc, cigma::array<float>& A);
void getDoubleArray(const char *loc, cigma::array<double>& A);
+ Cell::type getCellType(const char *loc);
+
boost::shared_ptr<NodeCoordinates> getNodeCoordinates(const char *loc);
boost::shared_ptr<ElementBlock> getElementBlock(const char *loc);
boost::shared_ptr<MeshPart> getMeshPart(const char *loc);
@@ -43,22 +46,22 @@
public:
std::string filename;
-
+ boost::shared_ptr<NcFile> ncfile;
};
#else
-/* If HAVE_VTK is not defined, we subclass NullReader */
+/* If HAVE_NETCDF is not defined, we subclass NullReader */
#include "io_null_reader.h"
-class cigma::ExodusII_Reader : public cigma::NullReader
+class cigma::ExodusReader : public cigma::NullReader
{
public:
- ExodusII_Reader();
- ExodusII_Reader();
+ ExodusReader();
+ ~ExodusReader();
ReaderType getReaderType() { return EXO_FILE_READER; }
};
#endif /* HAVE_VTK */
-#endif /* __CIGMA_EXODUSII_READER_H__ */
+#endif /* __CIGMA_EXODUS_READER_H__ */
Added: cs/cigma/trunk/src/io_netcdf.cpp
===================================================================
--- cs/cigma/trunk/src/io_netcdf.cpp (rev 0)
+++ cs/cigma/trunk/src/io_netcdf.cpp 2009-01-12 22:10:06 UTC (rev 13838)
@@ -0,0 +1,166 @@
+#include "io_netcdf.h"
+
+#ifdef HAVE_NETCDF
+
+#include <cassert>
+#include <stdexcept>
+#include <sstream>
+
+using namespace std;
+
+void nc_read_coordinates(NcFile& ncfile, valarray<double>& coordinates, int& num_vertices, int& num_dims)
+{
+ NcDim* num_nodes = ncfile.get_dim("num_nodes");
+ if (num_nodes == 0)
+ {
+ throw runtime_error("Could not read dimension 'num_nodes'");
+ }
+ num_vertices = num_nodes->size();
+
+ NcVar* coord = ncfile.get_var("coord");
+ if (coord == 0)
+ {
+ throw runtime_error("Could not get variable 'coord'");
+ }
+ if (coord->num_dims() != 2)
+ {
+ throw runtime_error("Number of dimensions of variable 'coord' must be 2.");
+ }
+
+ const int size = coord->num_vals();
+
+ NcDim* space_dim = coord->get_dim(0);
+ if (space_dim == 0)
+ {
+ throw runtime_error("Could not get dimensions of coordinates");
+ }
+ num_dims = space_dim->size();
+
+ assert(num_vertices * num_dims == size);
+
+ long* counts = coord->edges();
+ valarray<double> buffer(size);
+ bool ok = coord->get(&buffer[0], counts);
+ delete [] counts;
+ counts = 0;
+ if (!ok)
+ {
+ throw runtime_error("Could not get coordinate values");
+ }
+
+
+ coordinates.resize(size);
+
+ // transpose buffer
+ for (int i = 0; i < num_vertices; i++)
+ {
+ for (int j = 0; j < num_dims; j++)
+ {
+ coordinates[num_dims*i + j] = buffer[num_vertices*j + i];
+ }
+ }
+
+}
+
+void nc_read_connectivity(NcFile& ncfile, valarray<int>& cells, int& num_cells, int& num_corners)
+{
+ NcDim* num_elem = ncfile.get_dim("num_elem");
+ if (num_elem == 0)
+ {
+ throw runtime_error("Could not get dimension 'num_elem'");
+ }
+ num_cells = num_elem->size();
+
+ NcDim* num_el_blk = ncfile.get_dim("num_el_blk");
+ if (num_el_blk == 0)
+ {
+ throw runtime_error("Could not get dimension 'num_el_blk'");
+ }
+
+ NcVar* eb_prop1 = ncfile.get_var("eb_prop1");
+ if (eb_prop1 == 0)
+ {
+ throw runtime_error("Could not get variable 'eb_prop1'");
+ }
+
+ std::ostringstream varname;
+ varname << "num_nod_per_el" << 1;
+ NcDim* num_nod_per_el = ncfile.get_dim(varname.str().c_str());
+ if (num_nod_per_el == 0)
+ {
+ throw runtime_error("Could not get dimension 'num_node_per_el'");
+ }
+ num_corners = num_nod_per_el->size();
+ const int size = num_cells * num_corners;
+ cells.resize(size);
+
+ varname.str("");
+ varname << "num_el_in_blk" << 1;
+ NcDim* num_el_in_blk = ncfile.get_dim(varname.str().c_str());
+ if (num_el_in_blk == 0)
+ {
+ throw runtime_error("Could not get dimension 'num_el_in_blk'");
+ }
+ const int block_size = num_el_in_blk->size();
+
+ varname.str("");
+ varname << "connect" << 1;
+ NcVar* connect = ncfile.get_var(varname.str().c_str());
+ if (connect == 0)
+ {
+ throw runtime_error("Could not get variable 'connect'");
+ }
+ if (connect->num_dims() != 2)
+ {
+ throw runtime_error("Number of dimensions of variable 'connect' must be 2.");
+ }
+
+ //const int size = connect->num_vals();
+ assert(size == connect->num_vals());
+
+ long* counts = connect->edges();
+ bool ok = connect->get(&cells[0], counts);
+ delete [] counts; counts = 0;
+ if (!ok)
+ {
+ throw runtime_error("Could not get cell values");
+ }
+
+ // use zero-based indexing
+ cells -= 1;
+}
+
+void nc_read_elem_type(NcFile& ncfile, const char *varname, std::string& elem_type)
+{
+ NcVar* var = ncfile.get_var(varname);
+ if (var == 0)
+ {
+ elem_type = "";
+ }
+ else
+ {
+ NcAtt* att = var->get_att("elem_type");
+ if (att == 0)
+ {
+ elem_type = "";
+ }
+ else
+ {
+ NcValues* vals = att->values();
+ if (vals == 0)
+ {
+ elem_type = "";
+ }
+ else
+ {
+ const char *elt = vals->as_string(0);
+ elem_type = std::string(elt);
+ }
+ delete vals;
+ }
+ delete att;
+ }
+ return;
+}
+
+#endif /* HAVE_NETCDF */
Added: cs/cigma/trunk/src/io_netcdf.h
===================================================================
--- cs/cigma/trunk/src/io_netcdf.h (rev 0)
+++ cs/cigma/trunk/src/io_netcdf.h 2009-01-12 22:10:06 UTC (rev 13838)
@@ -0,0 +1,19 @@
+#ifndef CIGMA_IO_NETCDF_H
+#define CIGMA_IO_NETCDF_H
+
+#include "Common.h"
+
+#ifdef HAVE_NETCDF
+
+#include <valarray>
+#include <netcdfcpp.h>
+
+
+void nc_read_coordinates(NcFile& ncfile, std::valarray<double>& coordinates, int& num_vertices, int& num_dims);
+void nc_read_connectivity(NcFile& ncfile, std::valarray<int>& connectivity, int& num_cells, int& num_corners);
+void nc_read_elem_type(NcFile& ncfile, const char *varname, std::string& elem_type);
+
+
+#endif /* HAVE_NETCDF */
+
+#endif /* CIGMA_IO_NETCDF_H */
More information about the CIG-COMMITS
mailing list