XdmfHDF5Writer.cpp 29 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 <cmath>
28
#include <set>
29
#include <list>
30
#include "XdmfItem.hpp"
31
#include "XdmfArray.hpp"
32
#include "XdmfArrayType.hpp"
33
#include "XdmfError.hpp"
34
#include "XdmfHDF5Controller.hpp"
35
#include "XdmfHDF5Writer.hpp"
36

37 38
namespace {

39
  const static unsigned int DEFAULT_CHUNK_SIZE = 1000;
40 41 42

}

43 44 45
/**
 * PIMPL
 */
46
class XdmfHDF5Writer::XdmfHDF5WriterImpl  {
47 48 49 50

public:

  XdmfHDF5WriterImpl():
51
    mHDF5Handle(-1),
52
    mChunkSize(DEFAULT_CHUNK_SIZE),
53
    mOpenFile(""),
54
    mDepth(0)
55 56 57 58 59 60 61 62 63 64 65 66
  {
  };

  ~XdmfHDF5WriterImpl()
  {
    closeFile();
  };

  void
  closeFile()
  {
    if(mHDF5Handle >= 0) {
67
      /*herr_t status =*/H5Fclose(mHDF5Handle);
68 69
      mHDF5Handle = -1;
    }
70
    mOpenFile = "";
71 72
  };  

73 74
  int
  openFile(const std::string & filePath,
75 76
           const int fapl,
           const int mDataSetId)
77 78 79 80 81 82 83 84 85 86 87
  {
    if(mHDF5Handle >= 0) {
      // Perhaps we should throw a warning.
      closeFile();
    }
    // 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);
  
88
    int toReturn = 0;
89

90 91
    mOpenFile.assign(filePath);

92

93 94 95 96
    if(H5Fis_hdf5(filePath.c_str()) > 0) {
      mHDF5Handle = H5Fopen(filePath.c_str(), 
                            H5F_ACC_RDWR, 
                            fapl);
97 98 99 100 101 102 103 104 105
      if(mDataSetId == 0) {
        hsize_t numObjects;
        /*herr_t status = */H5Gget_num_objs(mHDF5Handle,
                                            &numObjects);
        toReturn = numObjects;
      }
      else {
        toReturn = mDataSetId;
      }
106 107
    }
    else {
108
      // This is where it currently fails
109 110 111 112 113 114 115 116 117
      mHDF5Handle = H5Fcreate(filePath.c_str(),
                              H5F_ACC_TRUNC,
                              H5P_DEFAULT,
                              fapl);
    }

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

118 119
    return toReturn;

120 121 122
  }

  hid_t mHDF5Handle;
123
  unsigned int mChunkSize;
124
  std::string mOpenFile;
125 126
  int mDepth;
  std::set<const XdmfItem *> mWrittenItems;
127 128

};
129

130
shared_ptr<XdmfHDF5Writer>
131 132
XdmfHDF5Writer::New(const std::string & filePath,
                    const bool clobberFile)
133
{
134 135 136
  if(clobberFile) {
    std::remove(filePath.c_str());
  }
137
  shared_ptr<XdmfHDF5Writer> p(new XdmfHDF5Writer(filePath));
138
  return p;
139 140
}

141
XdmfHDF5Writer::XdmfHDF5Writer(const std::string & filePath) :
142
  XdmfHeavyDataWriter(filePath, 1, 800),
143
  mImpl(new XdmfHDF5WriterImpl())
144 145 146 147 148
{
}

XdmfHDF5Writer::~XdmfHDF5Writer()
{
149
  delete mImpl;
150 151
}

152 153
shared_ptr<XdmfHeavyDataController>
XdmfHDF5Writer::createController(const std::string & hdf5FilePath,
154
                                     const std::string & dataSetPath,
155
                                     const shared_ptr<const XdmfArrayType> type,
156 157
                                     const std::vector<unsigned int> & start,
                                     const std::vector<unsigned int> & stride,
158 159
                                     const std::vector<unsigned int> & dimensions,
                                     const std::vector<unsigned int> & dataspaceDimensions)
160
{
161 162 163 164 165 166 167 168 169 170 171 172
  try {
    return XdmfHDF5Controller::New(hdf5FilePath,
                                   dataSetPath,
                                   type,
                                   start,
                                   stride,
                                   dimensions,
                                   dataspaceDimensions);
  }
  catch (XdmfError e) {
    throw e;
  }
173 174
}

175 176
unsigned int
XdmfHDF5Writer::getChunkSize() const
177
{
178
  return mImpl->mChunkSize;
179 180 181
}

int
182
XdmfHDF5Writer::getDataSetSize(std::string fileName, std::string dataSetName, const int fapl)
183
{
184 185 186 187 188 189 190
  hid_t handle = -1;
  H5E_auto_t old_func;
  void * old_client_data;
  H5Eget_auto(0, &old_func, &old_client_data);
  H5Eset_auto2(0, NULL, NULL);
  if (fileName !=  mImpl->mOpenFile) {
    // Save old error handler and turn off error handling for now
191

192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207
    if(H5Fis_hdf5(fileName.c_str()) > 0) {
      handle = H5Fopen(fileName.c_str(),
                       H5F_ACC_RDWR,
                       fapl);
    }
    else {
      // This is where it currently fails
      handle = H5Fcreate(fileName.c_str(),
                         H5F_ACC_TRUNC,
                         H5P_DEFAULT,
                         fapl);
    }
  }
  else {
    handle = mImpl->mHDF5Handle;
  }
208

209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225
  // Restore previous error handler
  H5Eset_auto2(0, old_func, old_client_data);

  hid_t checkset = H5Dopen(handle,
                           dataSetName.c_str(),
                           H5P_DEFAULT);
  hid_t checkspace = H5S_ALL;
  checkspace = H5Dget_space(checkset);
  hssize_t checksize = H5Sget_simple_extent_npoints(checkspace);
  herr_t status = H5Dclose(checkset);
  if(checkspace != H5S_ALL) {
    status = H5Sclose(checkspace);
  }
  if (handle != mImpl->mHDF5Handle) {
    H5Fclose(handle);
  }
  return checksize;
226 227
}

228 229 230 231 232 233 234 235 236
void 
XdmfHDF5Writer::closeFile()
{
  mImpl->closeFile();
}

void 
XdmfHDF5Writer::openFile()
{
237 238 239 240 241 242 243
  this->openFile(H5P_DEFAULT);
}

void
XdmfHDF5Writer::openFile(const int fapl)
{
  mDataSetId = mImpl->openFile(mFilePath,
244 245
                               fapl,
                               mDataSetId);
246 247
}

248 249 250 251 252 253
void
XdmfHDF5Writer::setChunkSize(const unsigned int chunkSize)
{
  mImpl->mChunkSize = chunkSize;
}

254 255
void
XdmfHDF5Writer::visit(XdmfArray & array,
256
                      const shared_ptr<XdmfBaseVisitor> visitor)
257
{
258 259 260
  mImpl->mDepth++;
  std::set<const XdmfItem *>::iterator checkWritten = mImpl->mWrittenItems.find(&array);
  if (checkWritten == mImpl->mWrittenItems.end() || array.getItemTag() == "DataItem") {
261
    // If it has children send the writer to them too.
262
    try {
263
      array.traverse(visitor);
264 265 266 267
    }
    catch (XdmfError e) {
      throw e;
    }
268
    if (array.isInitialized()) {
269
      // Only do this if the object has not already been written
270 271 272 273 274 275 276
      try {
        this->write(array, H5P_DEFAULT);
      }
      catch (XdmfError e) {
        throw e;
      }
      mImpl->mWrittenItems.insert(&array);
277
    }
278
  }
279
  // If the object has already been written, just end, it already has the data
280 281 282 283 284 285 286 287 288 289 290
  mImpl->mDepth--;
  if(mImpl->mDepth <= 0) {
    mImpl->mWrittenItems.clear();
  }
}

void
XdmfHDF5Writer::visit(XdmfItem & item,
                      const shared_ptr<XdmfBaseVisitor> visitor)
{
  mImpl->mDepth++;
291 292 293
  // This is similar to the algorithm for writing XPaths
  // shouldn't be a problem if XPaths are turned off because all this does is avoid writing an object twice
  // if it was written once then all instances of the object should have the controller
294 295 296
  std::set<const XdmfItem *>::iterator checkWritten = mImpl->mWrittenItems.find(&item);
  if (checkWritten == mImpl->mWrittenItems.end()) {
    mImpl->mWrittenItems.insert(&item);
297 298 299 300 301 302
    try {
      item.traverse(visitor);
    }
    catch (XdmfError e) {
      throw e;
    }
303 304 305 306 307
  }
  mImpl->mDepth--;
  if(mImpl->mDepth <= 0) {
    mImpl->mWrittenItems.clear();
  }
308 309
}

310

311 312 313
void
XdmfHDF5Writer::write(XdmfArray & array,
                      const int fapl)
314
{
315
  hid_t datatype = -1;
Kenneth Leiter's avatar
Kenneth Leiter committed
316
  bool closeDatatype = false;
317

318
  // Determining data type
319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346
  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;
    }
347 348
    else if(array.getArrayType() == XdmfArrayType::String()) {
      // Strings are a special case as they have mutable size
Kenneth Leiter's avatar
Kenneth Leiter committed
349 350 351 352
      datatype = H5Tcopy(H5T_C_S1);
      H5Tset_size(datatype, H5T_VARIABLE);
      closeDatatype = true;
    }
353
    else {
354 355 356 357 358 359 360 361
      try {
        XdmfError::message(XdmfError::FATAL,
                           "Array of unsupported type in "
                           "XdmfHDF5Writer::write");
      }
      catch (XdmfError e) {
        throw e;
      }
362 363 364
    }
  }

365 366
  herr_t status;

367 368
  if(datatype != -1) {
    std::string hdf5FilePath = mFilePath;
369

370 371 372 373 374 375 376 377 378 379 380 381 382
    size_t extIndex;
    std::string checkFileName;
    std::string checkFileExt;
    extIndex = hdf5FilePath.find_last_of(".");
    if (extIndex == std::string::npos) {
      checkFileName = hdf5FilePath;
      checkFileExt = "";
    }
    else {
      checkFileName = hdf5FilePath.substr(0, extIndex);
      checkFileExt = hdf5FilePath.substr(extIndex+1);
    }

383 384
    std::stringstream dataSetPath;

385 386
    std::vector<shared_ptr<XdmfHeavyDataController> > previousControllers;

387
    // Hold the controllers in order to base the new controllers on them
388
    for(unsigned int i = 0; i < array.getNumberHeavyDataControllers(); ++i) {
389
      previousControllers.push_back(array.getHeavyDataController(i));
390
    }
391

392 393
    // Remove controllers from the array
    // they will be replaced by the controllers created by this function.
394
    while(array.getNumberHeavyDataControllers() != 0) {
395
      array.removeHeavyDataController(array.getNumberHeavyDataControllers() -1);
396 397 398
    }


399 400

    if (previousControllers.size() == 0) {
401
      // Create a temporary controller if the array doesn't have one
402
      try {
403 404 405 406 407 408 409 410
        shared_ptr<XdmfHeavyDataController> tempDataController =
          this->createController(hdf5FilePath,
                                 "Data",
                                 array.getArrayType(),
                                 std::vector<unsigned int>(1, 0),
                                 std::vector<unsigned int>(1, 1),
                                 std::vector<unsigned int>(1, array.getSize()),
                                 std::vector<unsigned int>(1, array.getSize()));
411 412 413 414 415
        previousControllers.push_back(tempDataController);
      }
      catch (XdmfError e) {
        throw e;
      }
416
    }
417

418
    int controllerIndexOffset = 0;
419

420 421
    // It is assumed that the array will have at least one controller
    // if it didn't have one a temporary one was generated
422
    for(unsigned int i = 0; i < previousControllers.size(); ++i)
423
    {
424 425
      if (mMode == Append) {
        // Append only cares about the last controller, so add the rest back in
426
	for (; i < previousControllers.size() - 1; ++i) {
427 428 429 430
          array.insert(previousControllers[i]);
	}
      }

431 432 433 434 435 436
      std::list<std::string> filesWritten;
      std::list<shared_ptr<XdmfArray> > arraysWritten;
      std::list<std::vector<unsigned int> > startsWritten;
      std::list<std::vector<unsigned int> > stridesWritten;
      std::list<std::vector<unsigned int> > dimensionsWritten;
      std::list<std::vector<unsigned int> > dataSizesWritten;
437
      std::list<unsigned int> arrayOffsetsWritten;
438 439 440 441 442 443 444 445 446 447 448

      // Open a hdf5 dataset and write to it on disk.
      hsize_t size = array.getSize();

      // 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);

      bool startedloop = false;
449 450 451 452 453 454
      unsigned int origFileIndex = getFileIndex();
      while ((mMode == Hyperslab
              && i < previousControllers.size())
             || !startedloop) {
        // Hyperslab mode wants to assign all data using the current location
        // without writing until all data sets are determined
455 456 457 458 459

        startedloop = true;

        shared_ptr<XdmfHeavyDataController> heavyDataController =
          previousControllers[i];
460
        // Stats for the data currently stored in the array
461 462 463 464 465 466 467 468 469 470 471 472 473 474 475
    
        std::vector<unsigned int> dimensions;
        if (mMode != Hyperslab) {
          dimensions = array.getDimensions();
        }
        else {
	  dimensions = heavyDataController->getDimensions();
        }
        std::vector<unsigned int> dataspaceDimensions = dimensions;
        std::vector<unsigned int> start(dimensions.size(), 0);
        std::vector<unsigned int> stride(dimensions.size(), 1);

        if((mMode == Overwrite || mMode == Append || mMode == Hyperslab)
          && heavyDataController) {

476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498
          // If overwriting, appending, or writing to a hyperslab this
          // should be an hdf5 controller - cast it and pull info of
          // dataset we are writing to.
          shared_ptr<XdmfHDF5Controller> hdf5Controller = 
            shared_dynamic_cast<XdmfHDF5Controller>(heavyDataController);

          if(hdf5Controller) {
            // Write to the previous dataset
            dataSetPath.str(std::string());
            dataSetPath << hdf5Controller->getDataSetPath();
            hdf5FilePath = hdf5Controller->getFilePath();
            if(mMode == Hyperslab) {
              // Start, stride, and dataspace dimensions only matter
              // for hyperslab mode
              dataspaceDimensions = hdf5Controller->getDataspaceDimensions();
              start = hdf5Controller->getStart();
              stride = hdf5Controller->getStride();
            }
          }
          else {
            XdmfError::message(XdmfError::FATAL,
                               "Can only overwrite, append, or write a "
                               "hyperslab to a dataset of the same type");
499 500 501 502 503
          }
        }
        else {
          dataSetPath.str(std::string());
          dataSetPath << "Data" << mDataSetId;
504 505
        }

506 507
        // Check here for if the file would become
        // larger than the limit after the addition.
508
        // Then check subsequent files for the same limitation
509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526
        controllerSplitting(array,
                            fapl,
                            controllerIndexOffset,
                            heavyDataController,
                            checkFileName,
                            checkFileExt,
                            dataSetPath.str(),
                            dimensions,
                            dataspaceDimensions,
                            start,
                            stride,
                            filesWritten,
                            arraysWritten,
                            startsWritten,
                            stridesWritten,
                            dimensionsWritten,
                            dataSizesWritten,
                            arrayOffsetsWritten);
527

528 529 530
        if (mMode == Hyperslab)
        {
          i++;
531
          setFileIndex(origFileIndex);
532
        }
533

534 535
      }

536 537 538 539 540 541
      std::list<std::string>::iterator fileNameWalker = filesWritten.begin();
      std::list<shared_ptr<XdmfArray> >::iterator arrayWalker = arraysWritten.begin();
      std::list<std::vector<unsigned int> >::iterator startWalker = startsWritten.begin();
      std::list<std::vector<unsigned int> >::iterator strideWalker = stridesWritten.begin();
      std::list<std::vector<unsigned int> >::iterator dimensionWalker = dimensionsWritten.begin();
      std::list<std::vector<unsigned int> >::iterator dataSizeWalker = dataSizesWritten.begin();
542
      std::list<unsigned int>::iterator arrayOffsetWalker = arrayOffsetsWritten.begin();
543

544
      // Loop based on the amount of blocks split from the array.
545
      for (unsigned int writeIndex = 0; writeIndex < arraysWritten.size(); ++writeIndex) {
546

547 548
	// This is the section where the data is written to hdf5
	// If you want to change the writer to write to a different data format, do it here
549

550 551 552 553 554 555
        std::string curFileName = *fileNameWalker;
        shared_ptr<XdmfArray> curArray = *arrayWalker;
        std::vector<unsigned int> curStart = *startWalker;
        std::vector<unsigned int> curStride = *strideWalker;
        std::vector<unsigned int> curDimensions = *dimensionWalker;
        std::vector<unsigned int> curDataSize = *dataSizeWalker;
556
        unsigned int curArrayOffset = *arrayOffsetWalker;
557

558

559
	bool closeFile = false;
560 561
        // This is meant to open files if it isn't already opened by the write prior
        // If it wasn't open prior to writing it will be closed after writing
562 563 564 565 566
        if (mImpl->mOpenFile.compare(curFileName) != 0) {
          if(mImpl->mHDF5Handle < 0) {
            closeFile = true;
          }
          mImpl->openFile(curFileName,
567
                          fapl, mDataSetId);
568
        }
569

570 571 572 573 574 575 576 577 578 579 580 581 582 583
	htri_t testingSet = H5Lexists(mImpl->mHDF5Handle,
                                      dataSetPath.str().c_str(),
                                      H5P_DEFAULT);

        hid_t dataset = 0;

        if (testingSet == 0) {
          dataset = -1;
        }
        else {
          dataset = H5Dopen(mImpl->mHDF5Handle,
                            dataSetPath.str().c_str(),
                            H5P_DEFAULT);
        }
584

585
        // If default mode find a new data set to write to (keep
586
        // incrementing dataSetId)
587
        if(dataset >=0 && mMode == Default) {
588 589 590 591
          while(true) {
            dataSetPath.str(std::string());
            dataSetPath << "Data" << ++mDataSetId;
            if(!H5Lexists(mImpl->mHDF5Handle,
592
                          dataSetPath.str().c_str(),
593 594 595 596 597 598 599 600 601 602 603 604 605 606 607
                          H5P_DEFAULT)) {
              dataset = H5Dopen(mImpl->mHDF5Handle,
                                dataSetPath.str().c_str(),
                                H5P_DEFAULT);
              break;
            }
          }
        }

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

        hid_t dataspace = H5S_ALL;
        hid_t memspace = H5S_ALL;

608 609
        std::vector<hsize_t> current_dims(curDataSize.begin(),
                                          curDataSize.end());
610

611 612
        if(dataset < 0) {
          // If the dataset doesn't contain anything
613

614
          std::vector<hsize_t> maximum_dims(curDimensions.size(), H5S_UNLIMITED);
615
          // Create a new dataspace
616
          dataspace = H5Screate_simple(current_dims.size(),
617 618 619
                                       &current_dims[0],
                                       &maximum_dims[0]);
          hid_t property = H5Pcreate(H5P_DATASET_CREATE);
620 621 622 623 624 625

          const hsize_t totalDimensionsSize =
            std::accumulate(current_dims.begin(),
                            current_dims.end(),
                            1,
                            std::multiplies<hsize_t>());
626
          // The Nth root of the chunk size divided by the dimensions added together
627 628
          const double factor =
            std::pow(((double)mImpl->mChunkSize / totalDimensionsSize),
629 630
                     1.0 / current_dims.size());
          // The end result is the amount of slots alloted per unit of dimension
631 632
          std::vector<hsize_t> chunk_size(current_dims.begin(),
                                          current_dims.end());
633 634
	  if (mImpl->mChunkSize > 0) {
            // The chunk size won't do anything unless it's positive
635 636 637 638 639 640 641 642
            for(std::vector<hsize_t>::iterator iter = chunk_size.begin();
                iter != chunk_size.end(); ++iter) {
              *iter = (hsize_t)(*iter * factor);
              if(*iter == 0) {
                *iter = 1;
              }
            }
          }
643 644

          status = H5Pset_chunk(property, current_dims.size(), &chunk_size[0]);
645
          // Use that dataspace to create a new dataset
646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663
          dataset = H5Dcreate(mImpl->mHDF5Handle,
                              dataSetPath.str().c_str(),
                              datatype,
                              dataspace,
                              H5P_DEFAULT,
                              property,
                              H5P_DEFAULT);
          status = H5Pclose(property);
        }

        if(mMode == Append) {
          // Need to resize dataset to fit new data

          // Get size of old dataset
          dataspace = H5Dget_space(dataset);
          hssize_t datasize = H5Sget_simple_extent_npoints(dataspace);
          status = H5Sclose(dataspace);

664
          // Reset the datasize if the file or set is different
665
          if (curFileName != previousControllers[i]->getFilePath()) {
666 667 668 669 670 671 672
            datasize = 0;
          }
          if (dataSetPath.str() != previousControllers[i]->getDataSetPath()) {
            datasize = 0;
          }

          // Resize to fit size of old and new data.
673
          hsize_t newSize = curArray->getSize() + datasize;
674
          status = H5Dset_extent(dataset, &newSize);
675
          
676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691
          // Select hyperslab to write to.
          memspace = H5Screate_simple(1, &size, NULL);
          dataspace = H5Dget_space(dataset);
          hsize_t dataStart = datasize;
          status = H5Sselect_hyperslab(dataspace,
                                       H5S_SELECT_SET,
                                       &dataStart,
                                       NULL,
                                       &size,
                                       NULL);
        }
        else if(mMode == Overwrite) {
          // Overwriting - dataset rank must remain the same (hdf5 constraint)
          dataspace = H5Dget_space(dataset);

          const unsigned int ndims = H5Sget_simple_extent_ndims(dataspace);
692 693
          if(ndims != current_dims.size()) {
            try {
694 695 696 697
            XdmfError::message(XdmfError::FATAL,                            \
                               "Data set rank different -- ndims != "
                               "current_dims.size() -- in "
                               "XdmfHDF5Writer::write");
698 699 700 701 702
            }
            catch (XdmfError e) {
              throw e;
            }
          }
703 704 705 706 707 708 709 710 711

          status = H5Dset_extent(dataset, &current_dims[0]);
          dataspace = H5Dget_space(dataset);
        }
        else if(mMode == Hyperslab) {
          // Hyperslab - dataset rank must remain the same (hdf5 constraint)
          dataspace = H5Dget_space(dataset);

          const unsigned int ndims = H5Sget_simple_extent_ndims(dataspace);
712 713 714 715 716 717 718 719 720 721 722
          if(ndims != current_dims.size()) {
            try {
              XdmfError::message(XdmfError::FATAL,                            \
                                 "Data set rank different -- ndims != "
                                 "current_dims.size() -- in "
                                 "XdmfHDF5Writer::write");
            }
            catch (XdmfError e) {
              throw e;
            }
          }
723 724 725
          status = H5Dset_extent(dataset, &current_dims[0]);
          dataspace = H5Dget_space(dataset);

726 727 728 729 730 731 732 733 734



          std::vector<hsize_t> count(curDimensions.begin(),
                                     curDimensions.end());
          std::vector<hsize_t> currStride(curStride.begin(),
                                          curStride.end());
          std::vector<hsize_t> currStart(curStart.begin(),
                                         curStart.end());
735 736 737 738 739 740 741 742 743 744 745 746

          memspace = H5Screate_simple(count.size(),
                                      &(count[0]),
                                      NULL);
          status = H5Sselect_hyperslab(dataspace,
                                       H5S_SELECT_SET,
                                       &currStart[0],
                                       &currStride[0],
                                       &count[0],
                                       NULL) ;

          if(status < 0) {
747 748 749 750 751 752 753 754
            try {
              XdmfError::message(XdmfError::FATAL,
                                 "H5Dset_extent returned failure in "
                                 "XdmfHDF5Writer::write -- status: " + status);
            }
            catch (XdmfError e) {
              throw e;
            }
755 756 757 758
          }
        }

        status = H5Dwrite(dataset,
759
                          datatype,
760
                          memspace,
761
                          dataspace,
762
                          H5P_DEFAULT,
763 764
                          curArray->getValuesInternal());

765
        if(status < 0) {
766 767 768 769 770 771 772 773
          try {
            XdmfError::message(XdmfError::FATAL,
                               "H5Dwrite returned failure in XdmfHDF5Writer::write "
                               "-- status: " + status);
          }
          catch (XdmfError e) {
            throw e;
          }
774
        }
775

776 777 778
        if(dataspace != H5S_ALL) {
          status = H5Sclose(dataspace);
        }
779

780 781 782
        if(memspace != H5S_ALL) {
          status = H5Sclose(memspace);
        }
783

784
        status = H5Dclose(dataset);
785

786
	// This is causing a lot of overhead
787 788 789
        if(closeFile) {
          mImpl->closeFile();
        }
790

791 792 793
        if(mMode == Default) {
          ++mDataSetId;
        }
794

795 796
        // Attach a new controller to the array
        shared_ptr<XdmfHDF5Controller> newDataController =
797 798
          shared_ptr<XdmfHDF5Controller>();
        //This generates an empty pointer
799 800 801

        unsigned int newSize;
        if(mMode == Append) {
802
          // Find data size
803
          mImpl->openFile(curFileName,
804
                          fapl, mDataSetId);
805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824
          hid_t checkset = H5Dopen(mImpl->mHDF5Handle,
                                   dataSetPath.str().c_str(),
                                   H5P_DEFAULT);
          hid_t checkspace = H5S_ALL;
          checkspace = H5Dget_space(checkset);
          newSize = H5Sget_simple_extent_npoints(checkspace);
          status = H5Dclose(checkset);
	  if(checkspace != H5S_ALL) {
	    status = H5Sclose(checkspace);
          }
 
          std::vector<unsigned int> insertStarts;
          insertStarts.push_back(0);
          std::vector<unsigned int> insertStrides;
          insertStrides.push_back(1);
          std::vector<unsigned int> insertDimensions;
          insertDimensions.push_back(newSize);
          std::vector<unsigned int> insertDataSpaceDimensions;
          insertDataSpaceDimensions.push_back(newSize);

825
          try {
826 827 828 829 830 831 832 833
            newDataController = 
              shared_dynamic_cast<XdmfHDF5Controller>(this->createController(curFileName,
                                                      dataSetPath.str(),
                                                      curArray->getArrayType(),
                                                      insertStarts,
                                                      insertStrides,
                                                      insertDimensions,
                                                      insertDataSpaceDimensions));
834 835 836 837
          }
          catch (XdmfError e) {
            throw e;
          }
838
        }
839

840 841
        if(!newDataController) {
          // If the controller wasn't generated by append
842 843
          try {
            newDataController =
844 845 846 847 848 849 850
              shared_dynamic_cast<XdmfHDF5Controller>(this->createController(curFileName,
                                                      dataSetPath.str(),
                                                      curArray->getArrayType(),
                                                      curStart,
                                                      curStride,
                                                      curDimensions,
                                                      curDataSize));
851 852 853 854
          }
          catch (XdmfError e) {
            throw e;
          }
855
        }
856

857 858
        newDataController->setArrayOffset(curArrayOffset);

859 860
        array.insert(newDataController);

861 862 863 864 865 866
        fileNameWalker++;
        arrayWalker++;
        startWalker++;
        strideWalker++;
        dimensionWalker++;
        dataSizeWalker++;
867
        arrayOffsetWalker++;
868 869


870
      }
871

872
    }
873

874 875
    if(closeDatatype) {
      status = H5Tclose(datatype);
876
    }
877 878 879 880

    if(mReleaseData) {
      array.release();
    }
881
  }
882
}