diff --git a/c_tools/mock/loaders/sdf_loader.cpp b/c_tools/mock/loaders/sdf_loader.cpp index f2b0c32..55cb279 100644 --- a/c_tools/mock/loaders/sdf_loader.cpp +++ b/c_tools/mock/loaders/sdf_loader.cpp @@ -166,7 +166,7 @@ public: } if (load_flags & NEED_VELOCITY) { - const char *name_template[3] = { "x", "y", "z" }; + const char *name_template[3] = { "vx", "vy", "vz" }; for (int q = 0; q < 3; q++) { d->Vel[q] = new float[numPartToLoad]; diff --git a/c_tools/zobov2/voz1b1/voz1b1.cpp b/c_tools/zobov2/voz1b1/voz1b1.cpp index 0bad2ce..d9ff28a 100644 --- a/c_tools/zobov2/voz1b1/voz1b1.cpp +++ b/c_tools/zobov2/voz1b1/voz1b1.cpp @@ -1,3 +1,4 @@ +#include #include #include #include @@ -59,10 +60,12 @@ void BoxData::checkParticle(double *xyz, bool& in_main, bool& in_buf) in_buf = in_buf && (fabs(xyz[d]) < totwidth2[d]); in_main = in_main && (fabs(xyz[d]) <= width2[d]); } + assert(!in_main || in_buf); } void BoxData::prepareBox(const PositionData& pdata, int *b_id, int *numdivs) { + double BF = std::numeric_limits::max(); guard_added = false; for (int i = 0; i < 3; i++) @@ -82,7 +85,7 @@ void BoxData::prepareBox(const PositionData& pdata, int *b_id, int *numdivs) exit(0); } g[i] = (bf / 2.)*(1. + sqrt(1 - 2.*s[i]*s[i]/(bf*bf))); - cout << format("s[%d] = %f, bf = %f, g[%d] = %f.") % s[i] % i % bf % g[i] % i << endl; + cout << format("s[%d] = %f, bf = %f, g[%d] = %f.") % i% s[i] % bf % i % g[i] << endl; c[i] = b_id[i] * width[i]; } @@ -99,18 +102,20 @@ void BoxData::prepareBox(const PositionData& pdata, int *b_id, int *numdivs) double xyz[3] = { pdata.xyz[0][i], pdata.xyz[1][i], pdata.xyz[2][i] }; bool is_it_in_buf, is_it_in_main; - checkParticle(xyz, is_it_in_buf, is_it_in_main); + checkParticle(xyz, is_it_in_main, is_it_in_buf); if (is_it_in_buf) nvpbuf++; if (is_it_in_main) nvp++; + assert(nvp <= nvpbuf); } nvpbuf += 6*(NGUARD+1)*(NGUARD+1); /* number of guard points */ + nvpall = nvpbuf; - parts = new coordT[3*nvpbuf]; - orig = new pid_t[nvpbuf]; + parts = new coordT[3*nvpall]; + orig = new pid_t[nvpall]; if (parts == 0) { @@ -123,12 +128,12 @@ void BoxData::prepareBox(const PositionData& pdata, int *b_id, int *numdivs) exit(0); } - nvp = 0; nvpall = 0; /* nvp = number of particles without buffer */ + nvp = 0; /* nvp = number of particles without buffer */ for (int j = 0; j < 3; j++) { - xyz_min[j] = -std::numeric_limits::max(); - xyz_max[j] = std::numeric_limits::max(); + xyz_min[j] = BF; + xyz_max[j] = -BF; } for (pid_t i = 0; i < pdata.np; i++) { @@ -138,6 +143,7 @@ void BoxData::prepareBox(const PositionData& pdata, int *b_id, int *numdivs) checkParticle(xyz, is_it_in_main, is_it_in_buf); if (is_it_in_main) { + assert(nvp < nvpall); for (int j = 0; j < 3; j++) { parts[3*nvp+j] = xyz[j]; @@ -160,6 +166,7 @@ void BoxData::prepareBox(const PositionData& pdata, int *b_id, int *numdivs) if (is_it_in_buf && !is_it_in_main) { + assert(nvpbuf < nvpall); for (int j = 0; j < 3; j++) { parts[3*nvpbuf+j] = xyz[j]; @@ -173,7 +180,6 @@ void BoxData::prepareBox(const PositionData& pdata, int *b_id, int *numdivs) } cout << format("nvpbuf = %d") % nvpbuf << endl; cout << format("x: %f,%f; y: %f,%f; z:%f,%f\n") % xyz_min[0] % xyz_max[0] % xyz_min[1] % xyz_max[1] % xyz_min[2] % xyz_max[2] << endl; - nvpall = nvpbuf; double predict = 1; for (int j = 0; j < 3; j++) @@ -188,6 +194,8 @@ void BoxData::addGuardPoints() int rot_si[3][3] = { {1, 0, 0}, {1,0,0}, {0,1,0} }; int rot_sj[3][3] = { {0, 1, 0}, {0,0,1}, {0,0,1} }; + pid_t nvpcur = nvpbuf; + for (int k = 0; k < 3; k++) { @@ -196,17 +204,19 @@ void BoxData::addGuardPoints() { for (int j=0; j <= NGUARD; j++) { + assert(nvpcur < nvpall); /* Bottom */ for (int a = 0; a < 3; a++) - parts[3*nvpall+a] = -width2[a] + (realT)i * s[a] * rot_si[k][a] + (realT)j * s[a] * rot_sj[k][a] - rot_g[k][a] * g[a]; + parts[3*nvpcur+a] = -width2[a] + (realT)i * s[a] * rot_si[k][a] + (realT)j * s[a] * rot_sj[k][a] - rot_g[k][a] * g[a]; - nvpall++; + nvpcur++; + assert(nvpcur < nvpall); /* Top */ for (int a = 0; a < 3; a++) - parts[3*nvpall+a] = (2*rot_g[k][a]-1)*width2[a] + (realT)i * s[a] * rot_si[k][a] + (realT)j * s[a] * rot_sj[k][a] + rot_g[k][a] * g[a]; + parts[3*nvpcur+a] = (2*rot_g[k][a]-1)*width2[a] + (realT)i * s[a] * rot_si[k][a] + (realT)j * s[a] * rot_sj[k][a] + rot_g[k][a] * g[a]; - nvpall++; + nvpcur++; } } } @@ -218,15 +228,15 @@ void BoxData::findBoundary() for (int j = 0; j < 3; j++) { - xyz_min[j] = -BF; - xyz_max[j] = BF; + 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]); + xyz_max[j] = std::max(xyz_max[j], parts[3*i+j]); } } } @@ -271,7 +281,6 @@ int main(int argc, char *argv[]) { char *suffix, *outDir, *posfile; PARTADJ *adjs; float *vols; - pid_t *orig; double border, boxsize[3]; int b[3], numdiv[3]; double totalvol; @@ -375,7 +384,7 @@ int main(int argc, char *argv[]) { /* Get the adjacencies back to their original values */ for (pid_t i=0; i