diff --git a/c_tools/zobov2/voz1b1/voz1b1.cpp b/c_tools/zobov2/voz1b1/voz1b1.cpp index 63348dd..0bad2ce 100644 --- a/c_tools/zobov2/voz1b1/voz1b1.cpp +++ b/c_tools/zobov2/voz1b1/voz1b1.cpp @@ -26,6 +26,7 @@ bool checkParameters(int *numdiv, int *b) return true; } + struct BoxData { float width[3], totwidth[3]; @@ -40,6 +41,8 @@ struct BoxData void prepareBox(const PositionData& pdata, int *b_id, int *numdivs); void checkParticle(double *xyz, bool& in_main, bool& in_buf); + void addGuardPoints(); + void findBoundary(); }; void BoxData::checkParticle(double *xyz, bool& in_main, bool& in_buf) @@ -189,9 +192,9 @@ void BoxData::addGuardPoints() { /* Add guard points */ - for (i=0; i<=NGUARD; i++) + for (int i=0; i <= NGUARD; i++) { - for (j=0; j<=NGUARD; j++) + for (int j=0; j <= NGUARD; j++) { /* Bottom */ for (int a = 0; a < 3; a++) @@ -209,6 +212,52 @@ void BoxData::addGuardPoints() } } +void BoxData::findBoundary() +{ + double BF = std::numeric_limits::max(); + + for (int j = 0; j < 3; j++) + { + xyz_min[j] = -BF; + xyz_max[j] = BF; + } + + for (pid_t i = nvpbuf; i < nvpall; i++) { + for (int j = 0; j < 3; j++) + { + xyz_min[j] = std::min(xyz_min[j], parts[3*i+j]); + xyz_max[j] = std::min(xyz_max[j], parts[3*i+j]); + } + } +} + +void saveTesselation(const string& outfile, PositionData& pdata, BoxData& boxdata, PARTADJ *adjs, float *vols) +{ + ofstream out(outfile.c_str()); + if (!out) + { + cout << format("Unable to open %s") % outfile << endl; + exit(0); + } + out.write((char *)&pdata.np, sizeof(pid_t)); + out.write((char *)&boxdata.nvp, sizeof(pid_t)); + cout << format("nvp = %d") % boxdata.nvp << endl; + + /* Tell us where the original particles were */ + out.write((char *)boxdata.orig, sizeof(pid_t)*boxdata.nvp); + /* Volumes*/ + out.write((char *)vols,sizeof(float)*boxdata.nvp); + /* Adjacencies */ + for (pid_t i = 0; i < boxdata.nvp; i++) + { + out.write((char*)&adjs[i].nadj, sizeof(pid_t)); + if (adjs[i].nadj > 0) + out.write((char *)adjs[i].adj, adjs[i].nadj*sizeof(pid_t)); + else + (cout << "0").flush(); + } + out.close(); +} int main(int argc, char *argv[]) { PositionData pdata; @@ -270,26 +319,18 @@ int main(int argc, char *argv[]) { else boxdata.bf = 0.1; - boxdata.prepareBox(pdata); + boxdata.prepareBox(pdata, b, numdiv); pdata.destroy(); boxdata.addGuardPoints(); - adjs = new PARTADJ[np]; + adjs = new PARTADJ[boxdata.nvpall]; if (adjs == 0) { cout << "Unable to allocate adjs" << endl; return 0; } - xmin = BF; xmax = -BF; ymin = BF; ymax = -BF; zmin = BF; zmax = -BF; - for (i=nvpbuf;i xmax) xmax = parts[3*i]; - if (parts[3*i+1] < ymin) ymin = parts[3*i+1]; - if (parts[3*i+1] > ymax) ymax = parts[3*i+1]; - if (parts[3*i+2] < zmin) zmin = parts[3*i+2]; - if (parts[3*i+2] > zmax) zmax = parts[3*i+2]; - } + boxdata.findBoundary(); cout << format("Added guard points to total %d points (should be %d)") % boxdata.nvpall % (boxdata.nvpbuf + 6*(NGUARD+1)*(NGUARD+1)) << endl; @@ -297,7 +338,7 @@ int main(int argc, char *argv[]) { /* Do tesselation*/ printf("File read. Tessellating ...\n"); fflush(stdout); - exitcode = delaunadj(boxdata.parts, boxdata.nvp, boxdata.nvpbuf, boxdata.nvpall, &adjs); + int exitcode = delaunadj(boxdata.parts, boxdata.nvp, boxdata.nvpbuf, boxdata.nvpall, &adjs); if (exitcode != 0) { printf("Error while tesselating. Stopping here."); fflush(stdout); @@ -306,9 +347,9 @@ int main(int argc, char *argv[]) { /* Now calculate volumes*/ printf("Now finding volumes ...\n"); fflush(stdout); - vols = new float[nvp]; + vols = new float[boxdata.nvp]; - for (pid_t i = 0; i < nvp; i++) + for (pid_t i = 0; i < boxdata.nvp; i++) { /* Just the original particles * Assign adjacency coordinate array*/ /* Volumes */ @@ -316,7 +357,7 @@ int main(int argc, char *argv[]) { { for (int d = 0; d < 3; d++) { - deladjs[3*j + d] = parts[3*adjs[i].adj[j]+d] - parts[3*i+d]; + deladjs[3*j + d] = boxdata.parts[3*adjs[i].adj[j]+d] - boxdata.parts[3*i+d]; if (deladjs[3*j+d] < -boxdata.width2[d]) deladjs[3*j+d] += boxdata.width[d]; @@ -326,51 +367,28 @@ int main(int argc, char *argv[]) { } exitcode = vorvol(deladjs, points, intpoints, adjs[i].nadj, &(vols[i])); - vols[i] *= np/V0; + vols[i] *= pdata.np/pdata.V0; if ((i % 1000) == 0) cout << format("%d: %d, %f") % i % adjs[i].nadj % vols[i] << endl; } /* Get the adjacencies back to their original values */ - for (pid_t i=0; i 0) - out.write((char *)adjs[i].adj, adjs[i].nadj*sizeof(pid_t)); - else - (cout << "0").flush(); - } - out.close(); + saveTesselation(outfile, pdata, boxdata, adjs, vols); delete[] adjs; return(0); diff --git a/c_tools/zobov2/voz1b1/voz_io.cpp b/c_tools/zobov2/voz1b1/voz_io.cpp index 6bf033f..712bd5b 100644 --- a/c_tools/zobov2/voz1b1/voz_io.cpp +++ b/c_tools/zobov2/voz1b1/voz_io.cpp @@ -1,3 +1,4 @@ +#include #include #include #include "voz_io.hpp" @@ -12,15 +13,15 @@ bool PositionData::readFrom(const string& fname) try { UnformattedRead f(fname.c_str()); - + f.beginCheckpoint(); np = f.readInt32(); - f.endCheckPoint(); + f.endCheckpoint(); for (int j = 0; j < 3; j++) { xyz[j] = new float[np]; - + f.beginCheckpoint(); for (int p = 0; p < np; p++) xyz[j][p] = f.readReal32(); @@ -39,7 +40,7 @@ bool PositionData::readFrom(const string& fname) { return false; } - + return true; } @@ -49,13 +50,13 @@ void PositionData::findExtrema() for (int i = 0; i < 3; i++) { - xyz_min[i] = BF; + xyz_min[i] = BF; xyz_max[i] = -BF; } for (int i = 0; i < 3; i++) { - for (pid_t p = 0; p < np; p++) + for (pid_t p = 0; p < np; p++) { xyz_min[p] = min(xyz_min[p], xyz[p][i]); xyz_max[p] = max(xyz_max[p], xyz[p][i]); diff --git a/c_tools/zobov2/voz1b1/voz_io.hpp b/c_tools/zobov2/voz1b1/voz_io.hpp index 01205bb..77a1e5c 100644 --- a/c_tools/zobov2/voz1b1/voz_io.hpp +++ b/c_tools/zobov2/voz1b1/voz_io.hpp @@ -8,7 +8,8 @@ struct PositionData float *xyz[3]; pid_t np; float xyz_min[3], xyz_max[3]; - + float V0; + void destroy() { for (int j = 0; j < 3; j++)