diff --git a/src/fortran.cpp b/src/fortran.cpp index 367bf40..e497ff1 100644 --- a/src/fortran.cpp +++ b/src/fortran.cpp @@ -71,6 +71,16 @@ UnformattedRead::~UnformattedRead() delete f; } +int64_t UnformattedRead::position() const +{ + return f->tellg(); +} + +void UnformattedRead::seek(int64_t pos) +{ + f->seekg(pos, istream::beg); +} + // Todo implement primitive description void UnformattedRead::setOrdering(Ordering o) { diff --git a/src/fortran.hpp b/src/fortran.hpp index 81414ad..48b4a1d 100644 --- a/src/fortran.hpp +++ b/src/fortran.hpp @@ -95,6 +95,11 @@ namespace CosmoTool void skip(int64_t off) throw (InvalidUnformattedAccess); + int64_t position() const; + void seek(int64_t pos); + + void readOrderedBuffer(void *buffer, int size) + throw (InvalidUnformattedAccess); protected: bool swapOrdering; CheckpointSize cSize; @@ -102,8 +107,6 @@ namespace CosmoTool uint64_t checkPointAccum; std::ifstream *f; - void readOrderedBuffer(void *buffer, int size) - throw (InvalidUnformattedAccess); }; class UnformattedWrite: public FortranTypes diff --git a/src/loadGadget.cpp b/src/loadGadget.cpp index 474cf9e..9f77bb4 100644 --- a/src/loadGadget.cpp +++ b/src/loadGadget.cpp @@ -105,6 +105,11 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, boost::function0 readToSingle; long float_size; + if (GadgetFormat > 2) { + cerr << "Unknown gadget format" << endl; + return 0; + } + if (id >= 0) { int k = snprintf(0, 0, "%s.%d", fname, id)+1; char *out_fname = new char[k]; @@ -124,6 +129,29 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, } + typedef std::map BlockMap; + BlockMap blockTable; + + if (GadgetFormat == 2) { + int64_t startBlock = 0; + char block[5]; + int32_t blocksize; + + try { + while (true) { + f->beginCheckpoint(); + f->readOrderedBuffer(block, 4); + block[4] = 0; + blocksize = f->readInt32(); + f->endCheckpoint(); + blockTable[block] = f->position(); + f->skip(blocksize); + } + } catch (EndOfFileException&) {} + + f->seek(startBlock); + } + data = new SimuData; if (data == 0) { delete f; @@ -131,19 +159,26 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, } ssize_t NumPart = 0, NumPartTotal = 0; +#define ENSURE(name) { \ + if (GadgetFormat == 2) { \ + BlockMap::iterator iter = blockTable.find(name); \ + if (iter == blockTable.end()) { \ + std::cerr << "GADGET2: Cannot find block named '" << name << "'" << endl; \ + if (data) delete data; \ + delete f; \ + return 0; \ + } \ + f->seek(iter->second); \ + } \ + } + try { + ENSURE("HEAD"); loadGadgetHeader(f, h, data, id); - if (GadgetFormat == 1) - velmul = sqrt(h.time); - else if (GadgetFormat == 2) - velmul = 1/(h.time); - else { - cerr << "unknown gadget format" << endl; - abort(); - } + velmul = sqrt(h.time); if (h.flag_doubleprecision) { //cout << "Gadget format with double precision" << endl; @@ -167,10 +202,11 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, return 0; } + if (loadflags & NEED_TYPE) { int p = 0; - + data->type = new int[data->NumPart]; for (int k = 0; k < 6; k++) for (int n = 0; n < h.npart[k]; n++,p++) @@ -189,6 +225,7 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, try { + ENSURE("POS "); f->beginCheckpoint(); if (f->getBlockSize() != NumPart*float_size*3) { // Check that single would work @@ -239,6 +276,7 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, try { + ENSURE("VEL "); f->beginCheckpoint(); for(int k = 0, p = 0; k < 6; k++) { for(int n = 0; n < h.npart[k]; n++) { @@ -271,6 +309,7 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, if (loadflags & NEED_GADGET_ID) { try { + ENSURE("ID "); f->beginCheckpoint(); data->Id = new long[data->NumPart]; if (data->Id == 0) @@ -314,8 +353,10 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, try { long l = 0; - if (do_load) + if (do_load) { + ENSURE("MASS"); f->beginCheckpoint(); + } data->Mass = new float[NumPart]; for (int k = 0; k < 6; k++) { @@ -360,6 +401,7 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, return data; } +#undef ENSURE void CosmoTool::writeGadget(const std::string& fname, SimuData *data, int GadgetFormat) @@ -419,8 +461,7 @@ void CosmoTool::writeGadget(const std::string& fname, SimuData *data, int Gadget f->endCheckpoint(); float velmul = 1.0; - if (GadgetFormat == 1) - velmul = sqrt(data->time); + velmul = sqrt(data->time); f->beginCheckpoint(); for(int n = 0; n < data->NumPart; n++) {