diff --git a/src/loadGadget.cpp b/src/loadGadget.cpp index 38aff37..ee30d80 100644 --- a/src/loadGadget.cpp +++ b/src/loadGadget.cpp @@ -10,61 +10,41 @@ using namespace CosmoTool; using namespace std; -PurePositionData *CosmoTool::loadGadgetPosition(const char *fname) -{ - PurePositionData *data; - int p, n; - UnformattedRead f(fname); - GadgetHeader h; - - data = new PurePositionData; - f.beginCheckpoint(); - for (int i = 0; i < 6; i++) - h.npart[i] = f.readInt32(); - for (int i = 0; i < 6; i++) - h.mass[i] = f.readReal64(); - h.time = f.readReal64(); - h.redshift = f.readReal64(); - h.flag_sfr = f.readInt32(); - h.flag_feedback = f.readInt32(); - for (int i = 0; i < 6; i++) - h.npartTotal[i] = f.readInt32(); - h.flag_cooling = f.readInt32(); - h.num_files = f.readInt32(); - data->BoxSize = h.BoxSize = f.readReal64(); - h.Omega0 = f.readReal64(); - h.OmegaLambda = f.readReal64(); - h.HubbleParam = f.readReal64(); - f.endCheckpoint(true); - - data->NumPart = 0; - for(int k=0; k<5; k++) - data->NumPart += h.npart[k]; - - data->pos = new FCoordinates[data->NumPart]; - - f.beginCheckpoint(); - for(int k = 0, p = 0; k < 5; k++) { - for(int n = 0; n < h.npart[k]; n++) { - data->pos[p][0] = f.readReal32(); - data->pos[p][1] = f.readReal32(); - data->pos[p][2] = f.readReal32(); - p++; - } - } - f.endCheckpoint(); - - // Skip velocities - f.skip((long)data->NumPart*3+2*4); - // Skip ids - return data; +void loadGadgetHeader(UnformattedRead *f, GadgetHeader& h, SimuData *data, int id) +{ + f->beginCheckpoint(); + for (int i = 0; i < 6; i++) + h.npart[i] = f->readInt32(); + for (int i = 0; i < 6; i++) + h.mass[i] = f->readReal64(); + data->time = h.time = f->readReal64(); + h.redshift = f->readReal64(); + h.flag_sfr = f->readInt32(); + h.flag_feedback = f->readInt32(); + for (int i = 0; i < 6; i++) + h.npartTotal[i] = f->readInt32(); + h.flag_cooling = f->readInt32(); + h.num_files = f->readInt32(); + data->BoxSize = h.BoxSize = f->readReal64(); + data->Omega_M = h.Omega0 = f->readReal64(); + data->Omega_Lambda = h.OmegaLambda = f->readReal64(); + data->Hubble = h.HubbleParam = f->readReal64(); + f->endCheckpoint(true); + + long NumPart = 0, NumPartTotal = 0; + for(int k=0; k<6; k++) + { + NumPart += h.npart[k]; + NumPartTotal += (id < 0) ? h.npart[k] : h.npartTotal[k]; + } + data->NumPart = NumPart; + data->TotalNumPart = NumPartTotal; } - - - -SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, int loadflags, int GadgetFormat) +SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, + int loadflags, int GadgetFormat, + SimuFilter filter) { SimuData *data; int p, n; @@ -98,35 +78,11 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, int loadflags, i } long NumPart = 0, NumPartTotal = 0; - + try { - f->beginCheckpoint(); - for (int i = 0; i < 6; i++) - h.npart[i] = f->readInt32(); - for (int i = 0; i < 6; i++) - h.mass[i] = f->readReal64(); - data->time = h.time = f->readReal64(); - h.redshift = f->readReal64(); - h.flag_sfr = f->readInt32(); - h.flag_feedback = f->readInt32(); - for (int i = 0; i < 6; i++) - h.npartTotal[i] = f->readInt32(); - h.flag_cooling = f->readInt32(); - h.num_files = f->readInt32(); - data->BoxSize = h.BoxSize = f->readReal64(); - data->Omega_M = h.Omega0 = f->readReal64(); - data->Omega_Lambda = h.OmegaLambda = f->readReal64(); - data->Hubble = h.HubbleParam = f->readReal64(); - f->endCheckpoint(true); - - for(int k=0; k<6; k++) - { - NumPart += h.npart[k]; - NumPartTotal += (id < 0) ? h.npart[k] : h.npartTotal[k]; - } - data->NumPart = NumPart; - data->TotalNumPart = NumPartTotal; + loadGadgetHeader(f, h, data, id); + if (GadgetFormat == 1) velmul = sqrt(h.time); else if (GadgetFormat == 2) @@ -135,6 +91,9 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, int loadflags, i cerr << "unknown gadget format" << endl; abort(); } + + NumPart = data->NumPart; + NumPartTotal = data->TotalNumPart; } catch (const InvalidUnformattedAccess& e) { diff --git a/src/loadGadget.hpp b/src/loadGadget.hpp index 69f0ea6..be49ad4 100644 --- a/src/loadGadget.hpp +++ b/src/loadGadget.hpp @@ -5,10 +5,9 @@ #include "loadSimu.hpp" namespace CosmoTool { - - PurePositionData *loadGadgetPosition(const char *fname); - SimuData *loadGadgetMulti(const char *fname, int id, int flags, int GadgetFormat = 1); + SimuData *loadGadgetMulti(const char *fname, int id, int flags, + int GadgetFormat = 1, SimuFilter filter = 0); // Only single snapshot supported void writeGadget(const char *fname, SimuData *data, int GadgetFormat = 1); diff --git a/src/loadSimu.hpp b/src/loadSimu.hpp index a7db941..8e4c507 100644 --- a/src/loadSimu.hpp +++ b/src/loadSimu.hpp @@ -1,6 +1,8 @@ #ifndef __COSMOTOOLBOX_HPP #define __COSMOTOOLBOX_HPP +#include +#include namespace CosmoTool { @@ -9,9 +11,24 @@ namespace CosmoTool static const int NEED_VELOCITY = 4; static const int NEED_TYPE = 8; + struct SimuParticle + { + float Pos[3]; + float Vel[3]; + int type; + int id; + + bool flag_vel, flag_type, flag_id; + }; + + typedef bool (*SimuFilter)(const SimuParticle& p); + class SimuData { public: + typedef void (*FreeFunction)(void *); + typedef std::map > AttributeMap; + float BoxSize; float time; float Hubble; @@ -25,6 +42,9 @@ namespace CosmoTool float *Pos[3]; float *Vel[3]; int *type; + + AttributeMap attributes; + public: SimuData() : Id(0),NumPart(0),type(0) { Pos[0]=Pos[1]=Pos[2]=0; Vel[0]=Vel[1]=Vel[2]=0; } ~SimuData() @@ -40,7 +60,37 @@ namespace CosmoTool delete[] type; if (Id) delete[] Id; + + for (AttributeMap::iterator i = attributes.begin(); + i != attributes.end(); + ++i) + { + if (i->second.second) + i->second.second(i->second.first); + } } + + template + T *as(const std::string& n) + { + AttributeMap::iterator i = attributes.find(n); + if (i == attributes.end()) + return 0; + + return reinterpret_cast(i->first); + } + + void new_attribute(const std::string& n, void *p, FreeFunction free_func) + { + AttributeMap::iterator i = attributes.find(n); + if (i != attributes.end()) + { + if (i->second.second) + i->second.second(i->second.first); + } + attributes[n] = std::make_pair(p, free_func); + } + }; };