XdmfHDF5Writer.cpp 10.6 KB
Newer Older
Kenneth Leiter's avatar
Kenneth Leiter committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
/*****************************************************************************/
/*                                    XDMF                                   */
/*                       eXtensible Data Model and Format                    */
/*                                                                           */
/*  Id : XdmfHDF5Writer.cpp                                                  */
/*                                                                           */
/*  Author:                                                                  */
/*     Kenneth Leiter                                                        */
/*     kenneth.leiter@arl.army.mil                                           */
/*     US Army Research Laboratory                                           */
/*     Aberdeen Proving Ground, MD                                           */
/*                                                                           */
/*     Copyright @ 2011 US Army Research Laboratory                          */
/*     All Rights Reserved                                                   */
/*     See Copyright.txt 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.                                                 */
/*                                                                           */
/*****************************************************************************/
23

24
#include <hdf5.h>
25
#include <sstream>
26
#include <cstdio>
27
#include "XdmfArray.hpp"
28
#include "XdmfArrayType.hpp"
29
#include "XdmfHDF5Controller.hpp"
30
#include "XdmfHDF5Writer.hpp"
31
#include "XdmfError.hpp"
32

33
shared_ptr<XdmfHDF5Writer>
34 35
XdmfHDF5Writer::New(const std::string & filePath,
                    const bool clobberFile)
36
{
37 38 39
  if(clobberFile) {
    std::remove(filePath.c_str());
  }
40
  shared_ptr<XdmfHDF5Writer> p(new XdmfHDF5Writer(filePath));
41
  return p;
42 43
}

44
XdmfHDF5Writer::XdmfHDF5Writer(const std::string & filePath) :
45
  XdmfHeavyDataWriter(filePath)
46 47 48 49 50
{
}

XdmfHDF5Writer::~XdmfHDF5Writer()
{
51 52
}

53
shared_ptr<XdmfHDF5Controller>
54 55
XdmfHDF5Writer::createHDF5Controller(const std::string & hdf5FilePath,
                                     const std::string & dataSetPath,
56
                                     const shared_ptr<const XdmfArrayType> type,
57 58 59
                                     const std::vector<unsigned int> & start,
                                     const std::vector<unsigned int> & stride,
                                     const std::vector<unsigned int> & count)
60
{
61 62 63 64 65 66
  return XdmfHDF5Controller::New(hdf5FilePath,
                                 dataSetPath,
                                 type,
                                 start,
                                 stride,
                                 count);
67 68
}

69 70
void
XdmfHDF5Writer::visit(XdmfArray & array,
71
                      const shared_ptr<XdmfBaseVisitor> visitor)
72
{
73
  this->write(array, H5P_DEFAULT);
74 75
}

76 77 78
void
XdmfHDF5Writer::write(XdmfArray & array,
                      const int fapl)
79
{
80

81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111
  hid_t datatype = -1;

  if(array.isInitialized()) {
    if(array.getArrayType() == XdmfArrayType::Int8()) {
      datatype = H5T_NATIVE_CHAR;
    }
    else if(array.getArrayType() == XdmfArrayType::Int16()) {
      datatype = H5T_NATIVE_SHORT;
    }
    else if(array.getArrayType() == XdmfArrayType::Int32()) {
      datatype = H5T_NATIVE_INT;
    }
    else if(array.getArrayType() == XdmfArrayType::Int64()) {
      datatype = H5T_NATIVE_LONG;
    }
    else if(array.getArrayType() == XdmfArrayType::Float32()) {
      datatype = H5T_NATIVE_FLOAT;
    }
    else if(array.getArrayType() == XdmfArrayType::Float64()) {
      datatype = H5T_NATIVE_DOUBLE;
    }
    else if(array.getArrayType() == XdmfArrayType::UInt8()) {
      datatype = H5T_NATIVE_UCHAR;
    }
    else if(array.getArrayType() == XdmfArrayType::UInt16()) {
      datatype = H5T_NATIVE_USHORT;
    }
    else if(array.getArrayType() == XdmfArrayType::UInt32()) {
      datatype = H5T_NATIVE_UINT;
    }
    else {
112 113 114
      XdmfError::message(XdmfError::FATAL,
                         "Array of unsupported type in "
                         "XdmfHDF5Writer::write");
115 116 117 118 119
    }
  }

  if(datatype != -1) {
    std::string hdf5FilePath = mFilePath;
120

121 122 123 124 125 126 127 128 129 130 131 132
    std::stringstream dataSetPath;

    if((mMode == Overwrite || mMode == Append)
       && array.getHeavyDataController()) {
      // Write to the previous dataset
      dataSetPath << array.getHeavyDataController()->getDataSetPath();
      hdf5FilePath = array.getHeavyDataController()->getFilePath();
    }
    else {
      dataSetPath << "Data" << mDataSetId;
    }

133 134
    const std::vector<unsigned int> & dimensions = array.getDimensions();

135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154
    // Open a hdf5 dataset and write to it on disk.
    herr_t status;
    hsize_t size = array.getSize();
    hid_t hdf5Handle;

    // Save old error handler and turn off error handling for now
    H5E_auto_t old_func;
    void * old_client_data;
    H5Eget_auto(0, &old_func, &old_client_data);
    H5Eset_auto2(0, NULL, NULL);

    if(H5Fis_hdf5(hdf5FilePath.c_str()) > 0) {
      hdf5Handle = H5Fopen(hdf5FilePath.c_str(), H5F_ACC_RDWR, fapl);
    }
    else {
      hdf5Handle = H5Fcreate(hdf5FilePath.c_str(),
                             H5F_ACC_TRUNC,
                             H5P_DEFAULT,
                             fapl);
    }
155
    
156 157 158 159
    hid_t dataset = H5Dopen(hdf5Handle,
                            dataSetPath.str().c_str(),
                            H5P_DEFAULT);

160 161 162 163 164 165 166 167 168 169
    // if default mode find a new data set to write to (keep
    // incrementing dataSetId)
    while(dataset >= 0 && mMode == Default) {
      dataSetPath.str(std::string());
      dataSetPath << "Data" << ++mDataSetId;
      dataset = H5Dopen(hdf5Handle,
                        dataSetPath.str().c_str(),
                        H5P_DEFAULT);
    }

170 171 172
    hid_t dataspace = H5S_ALL;
    hid_t memspace = H5S_ALL;

173
    std::vector<hsize_t> current_dims(dimensions.begin(), dimensions.end());
174

175
    if(dataset < 0) {
176
      std::vector<hsize_t> maximum_dims(dimensions.size(), H5S_UNLIMITED);
177 178
      memspace = H5Screate_simple(dimensions.size(),
                                  &current_dims[0],
179
                                  &maximum_dims[0]);
180
      hid_t property = H5Pcreate(H5P_DATASET_CREATE);
181 182
      std::vector<hsize_t> chunk_size(dimensions.size(), 1024);
      status = H5Pset_chunk(property, dimensions.size(), &chunk_size[0]);
183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215
      dataset = H5Dcreate(hdf5Handle,
                          dataSetPath.str().c_str(),
                          datatype,
                          memspace,
                          H5P_DEFAULT,
                          property,
                          H5P_DEFAULT);
      status = H5Pclose(property);
    }
    else {
      // Need to resize dataset to fit new data
      if(mMode == Append) {
        // Get size of old dataset
        dataspace = H5Dget_space(dataset);
        hssize_t datasize = H5Sget_simple_extent_npoints(dataspace);
        status = H5Sclose(dataspace);

        // Resize to fit size of old and new data.
        hsize_t newSize = size + datasize;
        status = H5Dset_extent(dataset, &newSize);

        // Select hyperslab to write to.
        memspace = H5Screate_simple(1, &size, NULL);
        dataspace = H5Dget_space(dataset);
        hsize_t start = datasize;
        status = H5Sselect_hyperslab(dataspace,
                                     H5S_SELECT_SET,
                                     &start,
                                     NULL,
                                     &size,
                                     NULL) ;
      }
      else {
216 217 218 219
        // Overwriting - dataset rank must remain the same (hdf5 constraint)
        hid_t dataspace = H5Dget_space(dataset);

        const int ndims = H5Sget_simple_extent_ndims(dataspace);
220
        if(ndims != current_dims.size())
221 222 223 224
          XdmfError::message(XdmfError::FATAL, \
                             "Data set rank different -- ndims != "
                             "current_dims.size() -- in "
                             "XdmfHDF5Writer::write");
225 226 227 228

        status = H5Dset_extent(dataset, &current_dims[0]);

        if(status < 0) {
229 230 231
          XdmfError::message(XdmfError::FATAL,
                            "H5Dset_extent returned failure in "
                             "XdmfHDF5Writer::write -- status: " + status);
232
        }
233 234 235 236 237 238 239 240
      }
    }
    status = H5Dwrite(dataset,
                      datatype,
                      memspace,
                      dataspace,
                      H5P_DEFAULT,
                      array.getValuesInternal());
241 242

    if(status < 0) {
243 244 245
      XdmfError::message(XdmfError::FATAL, 
                         "H5Dwrite returned failure in XdmfHDF5Writer::write "
                         "-- status: " + status);
246 247
    }

248 249 250 251 252 253 254 255 256 257 258 259
    if(dataspace != H5S_ALL) {
      status = H5Sclose(dataspace);
    }
    if(memspace != H5S_ALL) {
      status = H5Sclose(memspace);
    }
    status = H5Dclose(dataset);
    status = H5Fclose(hdf5Handle);

    // Restore previous error handler
    H5Eset_auto2(0, old_func, old_client_data);

260
    // Attach a new controller to the array
261 262
    shared_ptr<XdmfHDF5Controller> newDataController =
      shared_ptr<XdmfHDF5Controller>();
263

264 265
    unsigned int newSize = array.getSize();
    if(mMode == Append && array.getHeavyDataController()) {
266 267 268 269 270 271 272 273 274
      newSize = array.getSize() + array.getHeavyDataController()->getSize();
      newDataController =
        this->createHDF5Controller(hdf5FilePath,
                                   dataSetPath.str(),
                                   array.getArrayType(),
                                   std::vector<unsigned int>(1, 0),
                                   std::vector<unsigned int>(1, 1),
                                   std::vector<unsigned int>(1, newSize));
    }
275

276
    if(mMode == Default || !array.getHeavyDataController()) {
277
      ++mDataSetId;
278
    }
279

280 281 282 283 284
    if(!newDataController) {
      newDataController =
        this->createHDF5Controller(hdf5FilePath,
                                   dataSetPath.str(),
                                   array.getArrayType(),
285 286 287 288
                                   std::vector<unsigned int>(dimensions.size(),
                                                             0),
                                   std::vector<unsigned int>(dimensions.size(),
                                                             1),
289 290
                                   dimensions);
    }
291
    array.setHeavyDataController(newDataController);
292 293 294 295

    if(mReleaseData) {
      array.release();
    }
296
  }
297
}