[cig-commits] r11483 - cs/benchmark/cigma/trunk/src

luis at geodynamics.org luis at geodynamics.org
Wed Mar 19 12:00:24 PDT 2008


Author: luis
Date: 2008-03-19 12:00:23 -0700 (Wed, 19 Mar 2008)
New Revision: 11483

Modified:
   cs/benchmark/cigma/trunk/src/HdfReader.cpp
   cs/benchmark/cigma/trunk/src/HdfReader.h
   cs/benchmark/cigma/trunk/src/Reader.cpp
   cs/benchmark/cigma/trunk/src/Reader.h
   cs/benchmark/cigma/trunk/src/TextReader.cpp
   cs/benchmark/cigma/trunk/src/TextReader.h
   cs/benchmark/cigma/trunk/src/VtkReader.cpp
   cs/benchmark/cigma/trunk/src/VtkReader.h
Log:
More robust handling of XML VTK files. Also made various improvements
to the reader classes.


Modified: cs/benchmark/cigma/trunk/src/HdfReader.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/HdfReader.cpp	2008-03-19 19:00:20 UTC (rev 11482)
+++ cs/benchmark/cigma/trunk/src/HdfReader.cpp	2008-03-19 19:00:23 UTC (rev 11483)
@@ -1,127 +1,142 @@
-#include "HdfReader.h"
-#include "h5io.h"
 #include <cstdlib>
 #include <cassert>
 
+#include "HdfReader.h"
+#include "HdfFile.h"
+#include "HdfDataset.h"
+
+using namespace cigma;
+
 // ---------------------------------------------------------------------------
 
-cigma::HdfReader::HdfReader()
+HdfReader::HdfReader()
 {
-    file_id = -1;
 }
 
-cigma::HdfReader::~HdfReader()
+HdfReader::~HdfReader()
 {
 }
 
 // ---------------------------------------------------------------------------
 
 
-int cigma::HdfReader::
-open(std::string filename)
+int HdfReader::open(const char *filename)
 {
-    file_id = h5io_file_open(filename.c_str(), "r");
-    assert(file_id >= 0);
-    return 0; // XXX: change return value instead of using assert
-}
+    h5.open(filename, "r");
 
-
-void cigma::HdfReader::
-close()
-{
-    herr_t status = H5Fclose(file_id);
-    if (status < 0)
+    if (h5.file_id < 0)
     {
-        // XXX: emit warning?
+        return -1;
     }
+
+    return 0;
 }
 
+int HdfReader::close()
+{
+    return h5.close();
+}
+
 // ---------------------------------------------------------------------------
 
-void cigma::HdfReader::
-get_dataset(const char *loc, double **data, int *num, int *dim)
+int HdfReader::get_dataset(const char *loc, double **data, int *num, int *dim)
 {
     int rank;
     hid_t type_id;
+    H5T_class_t typeclass;
     hid_t dataset_id;
+
     herr_t status;
     int ierr = -1;
 
-    dataset_id = h5io_dset_open(file_id, loc, &type_id, &rank, NULL, NULL);
-    assert(H5Tget_class(type_id) == H5T_FLOAT);
-    assert(dataset_id >= 0);
-    assert((rank == 1) || (rank == 2));
+    dataset_id = HdfDataset::open(h5.file_id, loc, &type_id, &rank, NULL, NULL);
+
+    if (dataset_id < 0)
+    {
+        return -1;
+    }
+
+    typeclass = H5Tget_class(type_id);
+    if (typeclass != H5T_FLOAT)
+    {
+        status = H5Dclose(dataset_id);
+        return -2;
+    }
+
+    if ((rank != 1) && (rank != 2))
+    {
+        status = H5Dclose(dataset_id);
+        return -3;
+    }
+
     status = H5Dclose(dataset_id);
-    assert(status >= 0); // XXX: emit warning?
 
     if (rank == 2)
     {
-        ierr = h5io_dset_read2(file_id, loc, H5T_NATIVE_DOUBLE,
-                               (void **)data, num, dim);
+        ierr = HdfDataset::read2(h5.file_id, loc, H5T_NATIVE_DOUBLE,
+                                 (void **)data, num, dim);
+        if (ierr < 0)
+        {
+            status = H5Dclose(dataset_id);
+            return -4;
+        }
     }
     else if (rank == 1)
     {
-        ierr = h5io_dset_read1(file_id, loc, H5T_NATIVE_DOUBLE,
-                               (void **)data, num);
+        ierr = HdfDataset::read1(h5.file_id, loc, H5T_NATIVE_DOUBLE,
+                                 (void **)data, num);
+        
+        if (ierr < 0)
+        {
+            status = H5Dclose(dataset_id);
+            return -4;
+        }
+
         *dim = 1;
     }
-    assert(ierr >= 0);
+    
+    return 0;
 }
 
-void cigma::HdfReader::
-get_coordinates(const char *loc, double **coordinates, int *nno, int *nsd)
+int HdfReader::get_coordinates(const char *loc, double **coordinates, int *nno, int *nsd)
 {
-    get_dataset(loc, coordinates, nno, nsd);
+    return this->get_dataset(loc, coordinates, nno, nsd);
 }
 
-void cigma::HdfReader::
-get_connectivity(const char *loc, int **connectivity, int *nel, int *ndofs)
+int HdfReader::get_connectivity(const char *loc, int **connectivity, int *nel, int *ndofs)
 {
     int rank;
     hid_t type_id;
+    H5T_class_t typeclass;
     hid_t connect_id;
     herr_t status;
+    int ierr;
 
-    connect_id = h5io_dset_open(file_id, loc, &type_id, &rank, NULL, NULL);
-    assert(H5Tget_class(type_id) == H5T_INTEGER);
-    assert(connect_id >= 0);
-    assert(rank == 2);
-    status = H5Dclose(connect_id);
-    if (status < 0)
+
+    connect_id = HdfDataset::open(h5.file_id, loc, &type_id, &rank, NULL, NULL);
+    if (connect_id < 0)
     {
-        // XXX: emit warning
+        return -1;
     }
 
-    int ierr;
-    ierr = h5io_dset_read2(file_id, loc, H5T_NATIVE_INT,
-                           (void **)connectivity, nel, ndofs);
-}
+    typeclass = H5Tget_class(type_id);
+    if (typeclass != H5T_INTEGER)
+    {
+        status = H5Dclose(connect_id);
+        return -2;
+    }
 
+    if (rank != 2)
+    {
+        status = H5Dclose(connect_id);
+        return -3;
+    }
 
-// ---------------------------------------------------------------------------
+    status = H5Dclose(connect_id);
 
-void cigma::HdfReader::
-get_mesh(MeshPart *meshPart, const char *loc)
-{
-    // XXX: read coords_loc & connect_loc from metadata
+    ierr = HdfDataset::read2(h5.file_id, loc, H5T_NATIVE_INT,
+                             (void **)connectivity, nel, ndofs);
 }
 
-void cigma::HdfReader::
-get_mesh(MeshPart *meshPart, const char *coords_loc, const char *connect_loc)
-{
-    assert(meshPart != 0);
-    get_coordinates(coords_loc, &(meshPart->coords), &(meshPart->nno), &(meshPart->nsd));
-    get_connectivity(connect_loc, &(meshPart->connect), &(meshPart->nel), &(meshPart->nel));
-}
 
-void cigma::HdfReader::
-get_dofs(DofHandler *dofHandler, const char *loc)
-{
-    assert(dofHandler != 0);
-    assert(dofHandler->meshPart != 0);
-    get_dataset(loc, &(dofHandler->dofs), &(dofHandler->nno), &(dofHandler->ndim));
-}
-
-
-
 // ---------------------------------------------------------------------------

Modified: cs/benchmark/cigma/trunk/src/HdfReader.h
===================================================================
--- cs/benchmark/cigma/trunk/src/HdfReader.h	2008-03-19 19:00:20 UTC (rev 11482)
+++ cs/benchmark/cigma/trunk/src/HdfReader.h	2008-03-19 19:00:23 UTC (rev 11483)
@@ -1,16 +1,11 @@
 #ifndef __HDF_READER_H__
 #define __HDF_READER_H__
 
-#include <string>
 #include "hdf5.h"
 #include "Reader.h"
+#include "HdfFile.h"
 
-#include "QuadraturePoints.h"
-#include "MeshPart.h"
-#include "DofHandler.h"
-#include "FE_Field.h"
 
-
 namespace cigma
 {
     class HdfReader;
@@ -27,26 +22,16 @@
     ReaderType getType() { return HDF_READER; }
 
 public:
-    int open(std::string filename);
-    void close();
+    int open(const char *filename);
+    int close();
 
 public:
-    void get_dataset(const char *loc, double **data, int *num, int *dim);
-    void get_coordinates(const char *loc, double **coordinates, int *nno, int *nsd);
-    void get_connectivity(const char *loc, int **connectivity, int *nel, int *ndofs);
+    int get_dataset(const char *loc, double **data, int *num, int *dim);
+    int get_coordinates(const char *loc, double **coordinates, int *nno, int *nsd);
+    int get_connectivity(const char *loc, int **connectivity, int *nel, int *ndofs);
 
 public:
-    void get_quadrature(QuadraturePoints *quadrature, const char *loc);
-    void get_quadrature(QuadraturePoints *quadrature, const char *points_loc, const char *weights_loc);
-
-public:
-    void get_field(FE_Field *field, const char *loc);
-    void get_mesh(MeshPart *meshPart, const char *coords_loc, const char *connect_loc);
-    void get_mesh(MeshPart *meshPart, const char *loc);
-    void get_dofs(DofHandler *dofHandler, const char *loc);
-
-public:
-    hid_t file_id;
+    HdfFile h5;
 };
 
 

Modified: cs/benchmark/cigma/trunk/src/Reader.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/Reader.cpp	2008-03-19 19:00:20 UTC (rev 11482)
+++ cs/benchmark/cigma/trunk/src/Reader.cpp	2008-03-19 19:00:23 UTC (rev 11483)
@@ -1,48 +1,38 @@
-#include <cassert>
+#include <string>
 #include <cstdlib>
+#include <cassert>
+
 #include "Reader.h"
+#include "NullReader.h"
 #include "HdfReader.h"
 #include "TextReader.h"
 #include "VtkReader.h"
-#include "VtkXmlReader.h"
 
 using namespace cigma;
 
 
 // ---------------------------------------------------------------------------
 
-void new_reader(cigma::Reader **reader, std::string ext)
+cigma::Reader* NewReader(const char *fileext)
 {
+    std::string ext = fileext;
+
     if (ext == ".h5")
     {
-        *reader = new HdfReader();
-        return;
+        return new HdfReader();
     }
 
     if (ext == ".txt")
     {
-        *reader = new TextReader();
-        return;
+        return new TextReader();
     }
 
-    if (ext == ".vtk")
+    if ((ext == ".vtk") || (ext == ".vts") || (ext == ".vtu"))
     {
-        // XXX: instantiate legacy vtk reader
-        *reader = new VtkReader();
-        return;
+        return new VtkReader();
     }
 
-    if ((ext == ".vts") || (ext == ".vtu"))
-    {
-        // XXX instantiate xml vtk reader
-        *reader = new VtkXmlReader();
-        if (ext == ".vtu")
-        {
-            // XXX: unstructured .vtu file reader still needed...so just fail for now
-            assert(false);
-        }
-        return;
-    }
+    return new NullReader(fileext);
 }
 
 // ---------------------------------------------------------------------------

Modified: cs/benchmark/cigma/trunk/src/Reader.h
===================================================================
--- cs/benchmark/cigma/trunk/src/Reader.h	2008-03-19 19:00:20 UTC (rev 11482)
+++ cs/benchmark/cigma/trunk/src/Reader.h	2008-03-19 19:00:23 UTC (rev 11483)
@@ -1,8 +1,6 @@
 #ifndef __READER_H__
 #define __READER_H__
 
-#include <string>
-
 namespace cigma
 {
     class Reader;
@@ -15,27 +13,27 @@
         NULL_READER,
         HDF_READER,
         VTK_READER,
-        TXT_READER
+        TEXT_READER
     } ReaderType;
 
+
 public:
     Reader();
     ~Reader();
 
 public:
     virtual ReaderType getType() = 0;
-    virtual int open(std::string filename) = 0;     // XXX: change std::string -> const char *
-    virtual void close() = 0;
 
 public:
-    // XXX: return error codes in the following functions
-    virtual void get_dataset(const char *loc, double **data, int *num, int *dim) = 0;
-    virtual void get_coordinates(const char *loc, double **coordinates, int *nno, int *nsd) = 0;
-    virtual void get_connectivity(const char *loc, int **connectivity, int *nel, int *ndofs) = 0;
+    virtual int open(const char *filename) = 0;
+    virtual int close() = 0;
+
+public:
+    virtual int get_dataset(const char *loc, double **data, int *num, int *dim) = 0;
+    virtual int get_coordinates(const char *loc, double **coordinates, int *nno, int *nsd) = 0;
+    virtual int get_connectivity(const char *loc, int **connectivity, int *nel, int *ndofs) = 0;
 };
 
+cigma::Reader* NewReader(const char *ext);
 
-void new_reader(cigma::Reader **reader, std::string ext);
-
-
 #endif

Modified: cs/benchmark/cigma/trunk/src/TextReader.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/TextReader.cpp	2008-03-19 19:00:20 UTC (rev 11482)
+++ cs/benchmark/cigma/trunk/src/TextReader.cpp	2008-03-19 19:00:23 UTC (rev 11483)
@@ -3,27 +3,34 @@
 #include <cstdlib>
 #include <cstring>
 
+
+using namespace cigma;
+
+
 // ---------------------------------------------------------------------------
-cigma::TextReader::TextReader()
+
+TextReader::TextReader()
 {
     fp = NULL;
 }
 
-cigma::TextReader::~TextReader()
+TextReader::~TextReader()
 {
     close();
 }
 
+
 // ---------------------------------------------------------------------------
-int cigma::TextReader::open(std::string filename)
+
+int TextReader::open(const char *filename)
 {
-    fp = NULL;
-
-    if (filename != "")
+    if (filename == NULL)
     {
-        fp = fopen(filename.c_str(), "r");
+        return -1;
     }
 
+    fp = fopen(filename, "r");
+
     // check for failure
     if (fp == NULL)
     {
@@ -33,12 +40,14 @@
     return 0;
 }
 
-void cigma::TextReader::close()
+int TextReader::close()
 {
     if (fp != NULL)
     {
         fclose(fp);
+        fp = NULL;
     }
+    return 0;
 }
 
 
@@ -122,32 +131,42 @@
     return true;
 }
 
-
-// ---------------------------------------------------------------------------
-
-void cigma::TextReader::get_dataset(const char *loc, double **data, int *num, int *dim)
+int TextReader::get_dataset(const char *loc, double **data, int *num, int *dim)
 {
     FILE *loc_fp = get_fp(loc, fp);
-    assert(loc_fp != NULL);
-    read_dmat(loc_fp, data, num, dim);
+    if (loc_fp != NULL)
+    {
+        read_dmat(loc_fp, data, num, dim);
+        return 0;
+    }
+    return -1;
 }
 
-void cigma::TextReader::get_connectivity(const char *loc, int **connectivity, int *nel, int *ndofs)
+int TextReader::get_connectivity(const char *loc, int **connectivity, int *nel, int *ndofs)
 {
     // XXX: add support for sections (c.f. gmsh format) so we can scan a single file
     // for now, interpret loc argument to be an entirely new file
     FILE *loc_fp = get_fp(loc, fp);
     assert(loc_fp != NULL);
-    read_imat(loc_fp, connectivity, nel, ndofs);
+    if (loc_fp != NULL)
+    {
+        read_imat(loc_fp, connectivity, nel, ndofs);
+        return 0;
+    }
+    return -1;
 }
 
-void cigma::TextReader::get_coordinates(const char *loc, double **coordinates, int *nno, int *nsd)
+int TextReader::get_coordinates(const char *loc, double **coordinates, int *nno, int *nsd)
 {
     // XXX: add support for sections (c.f. gmsh format) so we can scan a single file
     // for now, interpret loc argument to be an entirely new file
     FILE *loc_fp = get_fp(loc, fp);
-    assert(loc_fp != NULL);
-    read_dmat(loc_fp, coordinates, nno, nsd);
+    if (loc_fp != NULL)
+    {
+        read_dmat(loc_fp, coordinates, nno, nsd);
+        return 0;
+    }
+    return -1;
 }
 
 

Modified: cs/benchmark/cigma/trunk/src/TextReader.h
===================================================================
--- cs/benchmark/cigma/trunk/src/TextReader.h	2008-03-19 19:00:20 UTC (rev 11482)
+++ cs/benchmark/cigma/trunk/src/TextReader.h	2008-03-19 19:00:23 UTC (rev 11483)
@@ -1,16 +1,15 @@
-#ifndef __TEXTREADER_H__
-#define __TEXTREADER_H__
+#ifndef __TEXT_READER_H__
+#define __TEXT_READER_H__
 
 #include <cstdio>
-#include <string>
-
 #include "Reader.h"
 
 namespace cigma
 {
     class TextReader;
-};
+}
 
+
 class cigma::TextReader : public Reader
 {
 public:
@@ -18,16 +17,16 @@
     ~TextReader();
 
 public:
-    ReaderType getType() { return TXT_READER; }
+    ReaderType getType() { return TEXT_READER; }
 
 public:
-    int open(std::string filename);
-    void close();
+    int open(const char *filename);
+    int close();
 
 public:
-    void get_dataset(const char *loc, double **data, int *num, int *dim);
-    void get_coordinates(const char *loc, double **coordinates, int *nno, int *nsd);
-    void get_connectivity(const char *loc, int **connectivity, int *nel, int *ndofs);
+    int get_dataset(const char *loc, double **data, int *num, int *dim);
+    int get_coordinates(const char *loc, double **coordinates, int *nno, int *nsd);
+    int get_connectivity(const char *loc, int **connectivity, int *nel, int *ndofs);
 
 public:
     FILE *fp;   // default file pointer

Modified: cs/benchmark/cigma/trunk/src/VtkReader.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/VtkReader.cpp	2008-03-19 19:00:20 UTC (rev 11482)
+++ cs/benchmark/cigma/trunk/src/VtkReader.cpp	2008-03-19 19:00:23 UTC (rev 11483)
@@ -1,101 +1,149 @@
-#include <iostream>
+#include <string>
 #include <cstdlib>
 #include <cassert>
 
+#include "StringUtils.h"
+
 #include "VtkReader.h"
-
 #include "vtkCellArray.h"
 #include "vtkPointData.h"
 #include "vtkDoubleArray.h"
 #include "vtkFloatArray.h"
 
-
+
+using namespace std;
+using namespace cigma;
+
 // ---------------------------------------------------------------------------
 
-cigma::VtkReader::VtkReader()
+VtkReader::VtkReader()
 {
+    dataset = 0;
+    pointset = 0;
+    ugrid = 0;
+    sgrid = 0;
+
     reader = 0;
-    grid = 0;
+    xml_ug_reader = 0;
+    xml_sg_reader = 0;
 }
 
-cigma::VtkReader::~VtkReader()
+VtkReader::~VtkReader()
 {
-    //XXX: avoid calling close() for now
-    //close();
+    close();
 }
 
-
-
 // ---------------------------------------------------------------------------
 
-int cigma::VtkReader::open(std::string filename)
+
+int VtkReader::open(const char *filename)
 {
-    int ierr;
+    cout << "Reading " << filename << endl;
 
-    /* XXX: throw exception if file doesn't exist */
-    reader = vtkUnstructuredGridReader::New();
-    reader->SetFileName(filename.c_str());
+    string filepath = filename;
+    string fileroot, ext;
+    path_splitext(filepath, fileroot, ext);
 
+    if (ext == ".vtk")
+    {
+        int outputType;
 
-    //* detect bad vtk file
-    ierr = reader->OpenVTKFile();
-    if (ierr == 0)
+        reader = vtkDataSetReader::New();
+        reader->SetFileName(filename);
+        reader->Update();
+
+        outputType = reader->ReadOutputType();
+
+        if (outputType == VTK_STRUCTURED_GRID)
+        {
+            sgrid = reader->GetStructuredGridOutput();
+            dataset = sgrid;
+            pointset = sgrid;
+        }
+        else if (outputType == VTK_UNSTRUCTURED_GRID)
+        {
+            ugrid = reader->GetUnstructuredGridOutput();
+            dataset = ugrid;
+            pointset = ugrid;
+        }
+        else
+        {
+            cerr << "Unsupported VTK file. "
+                 << "Use only structured & unstructured grids."
+                 << endl;
+            return -1;
+        }
+    }
+    else if (ext == ".vts")
     {
-        //reader->Delete(); // XXX
-        reader = 0;
-        return -1;
+        xml_sg_reader = vtkXMLStructuredGridReader::New();
+        xml_sg_reader->SetFileName(filename);
+        xml_sg_reader->Update();
+
+        dataset = xml_sg_reader->GetOutputAsDataSet();
+        sgrid = vtkStructuredGrid::SafeDownCast(dataset);
     }
+    else if (ext == ".vtu")
+    {
+        xml_ug_reader = vtkXMLUnstructuredGridReader::New();
+        xml_ug_reader->SetFileName(filename);
+        xml_ug_reader->Update();
 
-    ierr = reader->ReadHeader();
-    if (ierr == 0)
+        dataset = xml_ug_reader->GetOutputAsDataSet();
+        pointset = vtkPointSet::SafeDownCast(dataset);
+        ugrid = vtkUnstructuredGrid::SafeDownCast(dataset);
+    }
+    else
     {
-        //reader->Delete(); // XXX
-        reader = 0;
+        cerr << "Unsupported VTK file type. "
+             << "Use only .vtk, .vts, or .vtu formats."
+             << endl;
         return -1;
-    } // */
+    }
 
-
-    reader->Update();
-    //reader->PrintSelf(std::cout, 0);
-
-    grid = reader->GetOutput();
-    //grid->PrintSelf(std::cout, 4);
-
     return 0;
 }
 
-void cigma::VtkReader::close()
+int VtkReader::close()
 {
-    /* XXX: fix this!
-    if (grid != 0)
+    if (this->reader)
     {
-        grid->Delete();
-        grid = 0;
+        this->reader->Delete();
+        this->reader = 0;
     }
-    if (reader != 0)
+    
+    if (this->xml_sg_reader)
     {
-        reader->Delete();
-        reader = 0;
+        this->xml_sg_reader->Delete();
+        this->xml_sg_reader = 0;
     }
-    // */
+
+    if (this->xml_ug_reader)
+    {
+        this->xml_ug_reader->Delete();
+        this->xml_ug_reader = 0;
+    }
+
+    return 0;
 }
 
 
-
 // ---------------------------------------------------------------------------
 
-void cigma::VtkReader::
-get_dataset(const char *loc, double **data, int *num, int *dim)
+
+int VtkReader::get_dataset(const char *loc, double **data, int *num, int *dim)
 {
-    assert(grid != 0);
-    assert(loc != 0);
-
     int size = 0;
     int dims[2] = {0, 0};
     double *array = 0;
 
-    vtkPointData *pointData = grid->GetPointData();
+    if (dataset == 0)
+    {
+        return -1;
+    }
 
+    vtkPointData *pointData = dataset->GetPointData();
+
     if (pointData->HasArray(loc) == 1)
     {
         vtkDataArray *dataArray = pointData->GetArray(loc);
@@ -127,160 +175,214 @@
         }
         else
         {
-            assert(false); // XXX: throw exception instead
+            assert(false); // XXX: throw own exception
+            return -1;
         }
     }
+    else
+    {
+        return -1; // no array found at loc!
+    }
 
+
     *data = array;
     *num = dims[0];
     *dim = dims[1];
-}
 
-void cigma::VtkReader::
-get_coordinates(const char *loc, double **coordinates, int *nno, int *nsd)
-{
-    get_coordinates(coordinates, nno, nsd);
-}
+    return 0;
 
-void cigma::VtkReader::
-get_connectivity(const char *loc, int **connectivity, int *nel, int *ndofs)
-{
-    get_connectivity(connectivity, nel, ndofs);
 }
 
-
-// ---------------------------------------------------------------------------
-
-void cigma::VtkReader::
-get_coordinates(double **coordinates, int *nno, int *nsd)
+int VtkReader::get_coordinates(const char *loc, double **coordinates, int *nno, int *nsd)
 {
-    assert(grid != 0);
+    if (ugrid != 0)
+    {
+        vtkPoints *points = pointset->GetPoints();
+        //points->PrintSelf(cout, 0);
 
-    vtkPoints *points = grid->GetPoints();
-    //points->PrintSelf(std::cout, 0);
+        int dims[2];
+        dims[0] = points->GetNumberOfPoints();
+        dims[1] = 3;
 
-    int dims[2];
-    dims[0] = points->GetNumberOfPoints();
-    dims[1] = 3;
+        int size = dims[0] * dims[1];
+        int dataType = points->GetDataType();
+        double *coords = 0;
 
-    int size = dims[0] * dims[1];
-    int dataType = points->GetDataType();
-    double *coords = new double[size];
-
-    if (dataType == VTK_DOUBLE)
-    {
-        double *ptr = static_cast<double*>(points->GetVoidPointer(0));
-        for (int i = 0; i < size; i++)
+        if (dataType == VTK_DOUBLE)
         {
-            coords[i] = ptr[i];
+            double *ptr;
+            coords = new double[size];
+            ptr = static_cast<double*>(points->GetVoidPointer(0));
+            for (int i = 0; i < size; i++)
+            {
+                coords[i] = ptr[i];
+            }
         }
-    }
-    else if (dataType == VTK_FLOAT)
-    {
-        float *ptr = static_cast<float*>(points->GetVoidPointer(0));
-        for (int i = 0; i < size; i++)
+        else if (dataType == VTK_FLOAT)
         {
-            coords[i] = ptr[i];
+            float *ptr;
+            coords = new double[size];
+            ptr = static_cast<float*>(points->GetVoidPointer(0));
+            for (int i = 0; i < size; i++)
+            {
+                coords[i] = ptr[i];
+            }
         }
+        else
+        {
+            assert(false); // XXX: throw own exception
+            return -1;
+        }
+
+        *coordinates = coords;
+        *nno = dims[0];
+        *nsd = dims[1];
+
+        return 0;
     }
-    else
+    else if (sgrid != 0)
     {
-        assert(false); // XXX: throw exception instead
-    }
+        vtkPoints *points = sgrid->GetPoints();
+        //points->PrintSelf(cout, 0);
 
-    *coordinates = coords;
-    *nno = dims[0];
-    *nsd = dims[1];
-}
 
-void cigma::VtkReader::
-get_connectivity(int **connectivity, int *nel, int *ndofs)
-{
-    assert(grid != 0);
+        // XXX: need to know if we're loading 2D cells, so we can
+        // recalculate size & nsd....if 3D, we leave things alone
 
-    vtkCellArray *cellArray = grid->GetCells();
-    //cellArray->PrintSelf(std::cout, 0);
 
-    vtkIdType numCells = grid->GetNumberOfCells();
-    vtkIdType *cellArrayPtr = cellArray->GetPointer();
+        int dims[2];
+        dims[0] = points->GetNumberOfPoints();
+        //dims[1] = 3;
+        dims[1] = 2;    // XXX: for now, assume 2D grid of quadrangles
 
-    int dofs_per_elt = cellArrayPtr[0];
-    int offset = dofs_per_elt + 1;
-    int *connect = new int[numCells * dofs_per_elt];
-    
-    for (int e = 0; e < numCells; ++e)
-    {
-        for (int i = 1; i <= dofs_per_elt; ++i)
+        int size = dims[0] * dims[1];
+        double *coords = new double[size];
+
+
+        int dataType = points->GetDataType();
+
+        if (dataType == VTK_DOUBLE)
         {
-            connect[dofs_per_elt * e + (i-1)] = cellArrayPtr[offset*e + i];
+            double *ptr = static_cast<double*>(points->GetVoidPointer(0));
+            for (int i = 0; i < dims[0]; i++)
+            {
+                coords[2*i + 0] = ptr[3*i + 0];
+                coords[2*i + 1] = ptr[3*i + 1];
+            } /*
+            for (int i = 0; i < size; i++)
+            {
+                coords[i] = ptr[i];
+            } // */
         }
+        else if (dataType == VTK_FLOAT)
+        {
+            float *ptr = static_cast<float*>(points->GetVoidPointer(0));
+            for (int i = 0; i < dims[0]; i++)
+            {
+                coords[2*i + 0] = ptr[3*i + 0];
+                coords[2*i + 1] = ptr[3*i + 1];
+            } /*
+            for (int i = 0; i < size; i++)
+            {
+                coords[i] = ptr[i];
+            } // */
+        }
+        else
+        {
+            assert(false); // XXX: throw exception instead
+        }
+
+        *coordinates = coords;
+        *nno = dims[0];
+        *nsd = dims[1];
+
+        return 0;
     }
 
-    *connectivity = connect;
-    *nel = numCells;
-    *ndofs = dofs_per_elt;
+    assert(false); // XXX: throw own exception
+    return -1;
 }
 
-void cigma::VtkReader::
-get_vector_point_data(const char *name, double **vectors, int *num, int *dim)
+
+int VtkReader::get_connectivity(const char *loc, int **connectivity, int *nel, int *ndofs)
 {
-    assert(grid != 0);
+    if (ugrid != 0)
+    {
+        vtkIdType numCells;
+        int dofs_per_elt;
+        int *connect = 0;
 
-    vtkPointData *pointData = grid->GetPointData();
-    //pointData->PrintSelf(std::cout, 0);
+        vtkCellArray *cellArray = ugrid->GetCells();
+        //cellArray->PrintSelf(std::cout, 0);
 
-    vtkDataArray *dataArray;
-    if (name != 0) {
-        assert(pointData->HasArray(name) == 1);
-        dataArray = pointData->GetVectors(name);
-    } else {
-        dataArray = pointData->GetVectors();
-    }
-    //dataArray->PrintSelf(std::cout, 0);
+        numCells = ugrid->GetNumberOfCells();
 
-    int dataType = dataArray->GetDataType();
-    assert(dataType == VTK_DOUBLE);
-    double *ptr = static_cast<double*>(dataArray->GetVoidPointer(0));
+        vtkIdType *cellArrayPtr = cellArray->GetPointer();
 
-    // XXX: copy the data, or keep the reference?
-    // if dataType from the file is float, then we'd need to copy anyway
-    // perhaps we need a void pointer and a type
-    // new function signature would be
-    //  void get_vector_point_data(const char *name, void **vectors, int *num, int *dim, int *type)
+        dofs_per_elt = cellArrayPtr[0];
 
-    *vectors = ptr;
-    *num = dataArray->GetNumberOfTuples();
-    *dim = dataArray->GetNumberOfComponents();
-}
+        int offset = dofs_per_elt + 1;
 
-void cigma::VtkReader::
-get_scalar_point_data(const char *name, double **scalars, int *num, int *dim)
-{
-    assert(grid != 0);
+        connect = new int[numCells * dofs_per_elt];
 
-    vtkPointData *pointData = grid->GetPointData();
-    //pointData->PrintSelf(std::cout, 0);
+        for (int e = 0; e < numCells; ++e)
+        {
+            for (int i = 1; i <= dofs_per_elt; ++i)
+            {
+                connect[dofs_per_elt * e + (i-1)] = cellArrayPtr[offset*e + i];
+            }
+        }
 
-    vtkDataArray *dataArray;
-    if (name != 0) {
-        assert(pointData->HasArray(name) == 1);
-        dataArray = pointData->GetScalars(name);
-    } else {
-        dataArray = pointData->GetScalars();
+        *connectivity = connect;
+        *nel = numCells;
+        *ndofs = dofs_per_elt;
+
+        return 0;
+
     }
-    //dataArray->PrintSelf(std::cout, 0);
+    else if (sgrid != 0)
+    {
+        //
+        // XXX: for now, assume 2D grid of quadrangles
+        //
+        int nx, ny;
+        int ex, ey;
+        int *connect = 0;
+        int ndofs_per_cell;
+        int num_cells;
+        int extent[6];
 
+        num_cells = sgrid->GetNumberOfCells();
+        ndofs_per_cell = 4;
 
-    int dataType = dataArray->GetDataType();
-    assert(dataType == VTK_DOUBLE);
-    double *ptr = static_cast<double*>(dataArray->GetVoidPointer(0));
+        sgrid->GetExtent(extent);
 
-    // XXX: see comment in get_vector_point_data()
+        nx = extent[1] + 1;
+        ny = extent[3] + 1;
 
-    *scalars = ptr;
-    *num = dataArray->GetNumberOfTuples();
-    *dim = dataArray->GetNumberOfComponents();
+        ex = nx - 1;
+        ey = ny - 1;
+
+        connect = new int[ex * ey * ndofs_per_cell];
+        for (int j = 0; j < ny-1; j++)
+        {
+            for (int i = 0; i < nx-1; i++)
+            {
+                int e = i + j*ex;
+                connect[4*e + 0] =  i    +  j*nx;
+                connect[4*e + 1] = (i+1) +  j*nx;
+                connect[4*e + 2] = (i+1) + (j+1)*nx;
+                connect[4*e + 3] =  i    + (j+1)*nx;
+            }
+        }
+
+        *connectivity = connect;
+        *nel = num_cells;
+        *ndofs = ndofs_per_cell;
+    }
+
+    assert(false); // XXX: throw own exception
+    return -1;
 }
 
-
 // ---------------------------------------------------------------------------

Modified: cs/benchmark/cigma/trunk/src/VtkReader.h
===================================================================
--- cs/benchmark/cigma/trunk/src/VtkReader.h	2008-03-19 19:00:20 UTC (rev 11482)
+++ cs/benchmark/cigma/trunk/src/VtkReader.h	2008-03-19 19:00:23 UTC (rev 11483)
@@ -1,80 +1,55 @@
 #ifndef __VTK_READER_H__
 #define __VTK_READER_H__
 
-#include <string> 
-
 #include "Reader.h"
 
 #include "vtkDataSetReader.h"
-#include "vtkUnstructuredGridReader.h"
+#include "vtkXMLStructuredGridReader.h"
+#include "vtkXMLUnstructuredGridReader.h"
+
+#include "vtkDataSet.h"
+#include "vtkPointSet.h"
+#include "vtkStructuredGrid.h"
 #include "vtkUnstructuredGrid.h"
 
+
+
 namespace cigma
 {
     class VtkReader;
 }
 
-class cigma::VtkReader : public cigma::Reader
+class cigma::VtkReader : public Reader
 {
-public:
-    typedef enum {
-        VTK_FILE_NONE,
-        VTK_FILE_POLYDATA,
-        VTK_FILE_STRUCTURED_POINTS,
-        VTK_FILE_STRUCTURED_GRID,
-        VTK_FILE_UNSTRUCTURED_GRID,
-        VTK_FILE_RECTILINEAR_GRID
-    } VtkFileType;
 
-    ReaderType getType()
-    {
-        return VTK_READER;
-    }
-
-    VtkFileType getFileType()
-    {
-        if (reader != 0)
-        {
-            // XXX: are these mutually exclusive?
-            if (reader->IsFilePolyData())
-                return VTK_FILE_POLYDATA;
-            if (reader->IsFileStructuredPoints())
-                return VTK_FILE_STRUCTURED_POINTS;
-            if (reader->IsFileStructuredGrid())
-                return VTK_FILE_STRUCTURED_GRID;
-            if (reader->IsFileUnstructuredGrid())
-                return VTK_FILE_UNSTRUCTURED_GRID;
-            if (reader->IsFileRectilinearGrid())
-                return VTK_FILE_RECTILINEAR_GRID;
-        }
-        return VTK_FILE_NONE;
-    }
-
 public:
     VtkReader();
     ~VtkReader();
 
 public:
-    int open(std::string filename);
-    void close();
+    ReaderType getType() { return VTK_READER; }
 
 public:
-    void get_dataset(const char *loc, double **data, int *num, int *dim);
-    void get_coordinates(const char *loc, double **coordinates, int *nno, int *nsd);
-    void get_connectivity(const char *loc, int **connectivity, int *nel, int *ndofs);
+    int open(const char *filename);
+    int close();
 
 public:
-    void get_coordinates(double **coordinates, int *nno, int *nsd);
-    void get_connectivity(int **connectivity, int *nel, int *ndofs);
-    void get_vector_point_data(const char *name, double **vectors, int *num, int *dim);
-    void get_scalar_point_data(const char *name, double **scalars, int *num, int *dim);
+    int get_dataset(const char *loc, double **data, int *num, int *dim);
+    int get_coordinates(const char *loc, double **coordinates, int *nno, int *nsd);
+    int get_connectivity(const char *loc, int **connectivity, int *nel, int *ndofs);
 
 public:
-    vtkUnstructuredGridReader *reader;
-    //vtkDataSetReader *reader;
-    vtkUnstructuredGrid *grid;
+    vtkDataSet *dataset;
+    vtkPointSet *pointset;
+    vtkUnstructuredGrid *ugrid;
+    vtkStructuredGrid *sgrid;
+
+public:
+    vtkDataSetReader *reader;
+    vtkXMLStructuredGridReader *xml_sg_reader;
+    vtkXMLUnstructuredGridReader *xml_ug_reader;
+
 };
 
-// ---------------------------------------------------------------------------
 
 #endif



More information about the cig-commits mailing list