From 53483fd0d881b3ff121b22a95ceb0adf589149be Mon Sep 17 00:00:00 2001 From: "Paul M. Sutter" Date: Sat, 28 Jun 2025 12:19:33 -0400 Subject: [PATCH] begin branch to replace qhull with voro++ for tessellation --- c_source/CMakeLists.txt | 2 +- c_source/vorozobov/.CMakeLists.txt.swp | Bin 0 -> 12288 bytes c_source/vorozobov/CMakeLists.txt | 15 + c_source/vorozobov/findrtop.c | 81 +++++ c_source/vorozobov/readfiles.cpp | 121 +++++++ c_source/vorozobov/voz.h | 11 + c_source/vorozobov/voz1b1.cpp | 481 +++++++++++++++++++++++++ c_source/vorozobov/vozinit.cpp | 203 +++++++++++ c_source/vorozobov/voztie.cpp | 285 +++++++++++++++ c_source/vorozobov/vozutil.cpp | 250 +++++++++++++ external/external_build.cmake | 1 + 11 files changed, 1449 insertions(+), 1 deletion(-) create mode 100644 c_source/vorozobov/.CMakeLists.txt.swp create mode 100644 c_source/vorozobov/CMakeLists.txt create mode 100644 c_source/vorozobov/findrtop.c create mode 100644 c_source/vorozobov/readfiles.cpp create mode 100644 c_source/vorozobov/voz.h create mode 100644 c_source/vorozobov/voz1b1.cpp create mode 100644 c_source/vorozobov/vozinit.cpp create mode 100644 c_source/vorozobov/voztie.cpp create mode 100644 c_source/vorozobov/vozutil.cpp diff --git a/c_source/CMakeLists.txt b/c_source/CMakeLists.txt index b36c0a6..e4b75b7 100644 --- a/c_source/CMakeLists.txt +++ b/c_source/CMakeLists.txt @@ -44,4 +44,4 @@ add_genopt(computeAverageDistortion_SRCS computeAverageDistortion.ggo computeAve #add_executable(computeAverageDistortion ${computeAverageDistortion_SRCS}) #target_link_libraries(computeAverageDistortion ${ZOB_LIBS}) -subdirs(util prep pruning jozov2 zobov) +subdirs(util prep pruning jozov2 vorozobov) diff --git a/c_source/vorozobov/.CMakeLists.txt.swp b/c_source/vorozobov/.CMakeLists.txt.swp new file mode 100644 index 0000000000000000000000000000000000000000..4284854ce4cab97dad2fb4e7acf26f767c300611 GIT binary patch literal 12288 zcmeI2v2W8r6vnR=Bc&BgY=&fL6eUhz0l|{A5w(~zcfa#u`NGD~>p73>O-5~vvF~fI_CI{S zvHbN3W65Eb8zWLXbWz*pxyoui9<{YTs4E?DHTUs1vx4qe#v(}aP|3)Rf=uUWB;2u1 z^|221*xl&y1JRY4$?E3F%&mo+AOHfFAds8A+M^Zr@IkX-vr6~xId^XFUZMaV1pyEM z0T2KI5C8!X009uV?gUJ-#NN^(mWzd~7IEpXSacCU00ck)1V8`;KmY_l00ck)1V8`; zt|0*)Gxp;qW8Mmxr~m)6@Bd$P@rCrB)Fd@XKW{Pijbum}X@hi^^oepmlHQSQou?!; z5C8!X009sH0T2KI5C8!X_^$}avu%_auGG5U8a(UxJI+`iH^Zh4awLt7Mha)Mv*UaH zmhV09JLSpH+w25wuYXcuC3!x7nz%57l?+c`3ll{8Fp)|GIx$j@vh`xB?SVIJ`5mX` zxp|hl(_`7L97QURPtRIct4^(ZW=zS@sVRF3yIQTVqPUWygMda-o=Wkjc6rik`CCEP zYm;>z$ALH!QEqss#BBOg7xgGfoK*05j~r&USD~|~i;Yj0d{xe;#s6V_wpljLHNV)q p-R`ThZfln(FSd5OT}yq?8oWHA8yD_E@%m9XzYzOA*7s>Iegn#4E&%`l literal 0 HcmV?d00001 diff --git a/c_source/vorozobov/CMakeLists.txt b/c_source/vorozobov/CMakeLists.txt new file mode 100644 index 0000000..8578d2b --- /dev/null +++ b/c_source/vorozobov/CMakeLists.txt @@ -0,0 +1,15 @@ + +add_executable(voz1b1 voz1b1.cpp readfiles.cpp voz.h) +target_link_libraries(voz1b1 ${QHULL_LIBRARY} ${MATH_LIB}) +target_compile_options(voz1b1 PRIVATE -I/usr/local/include/voro++ -L/usr/local/lib -lvoro++) + +add_executable(vozinit vozinit.cpp readfiles.cpp) +target_link_libraries(vozinit ${MATH_LIB}) +#target_compile_options(vozinit PRIVATE -I/usr/local/include/voro++ -L/usr/local/lib -lvoro++) + +add_executable(voztie voztie.cpp readfiles.cpp) +target_link_libraries(voztie ${MATH_LIB}) +#target_compile_options(voztie PRIVATE -I/usr/local/include/voro++ -L/usr/local/lib -lvoro++) + +#install(TARGETS vozinit DESTINATION ${VIDE_BIN}) +install(TARGETS voz1b1 vozinit voztie DESTINATION ${VIDE_BIN}) diff --git a/c_source/vorozobov/findrtop.c b/c_source/vorozobov/findrtop.c new file mode 100644 index 0000000..8dbeeb3 --- /dev/null +++ b/c_source/vorozobov/findrtop.c @@ -0,0 +1,81 @@ +#include +#include + +/*------------------------------------------------------------------------------ + Find nb elements of Real array a having the largest value. + Returns index iord of these elements, ordered so iord[0] corresponds + to element a[iord[0]] having the largest value. + If nb > na, then the last nb-na elements of iord are undefined. + Elements of a that are equal are left in their original order. + + Courtesy Andrew Hamilton. +*/ +void findrtop(double *a, int na, int *iord, int nb) +{ +#undef ORDER +#define ORDER(ia, ja) a[ia] > a[ja] || (a[ia] == a[ja] && ia < ja) + + int i, ia, ib, it, n; + + n = (na <= nb)? na : nb; + if (n <= 0) return; + + /* heap first n elements, so smallest element is at top of heap */ + for (ib = (n >> 1); ib < n; ib++) { + iord[ib] = ib; + } + for (ia = (n >> 1) - 1; ia >= 0; ia--) { + i = ia; + for (ib = (i << 1) + 1; ib < n; ib = (i << 1) + 1) { + if (ib+1 < n) { + if (ORDER(iord[ib], iord[ib+1])) ib++; + } + if (ORDER(ia, iord[ib])) { + iord[i] = iord[ib]; + i = ib; + } else { + break; + } + } + iord[i] = ia; + } + + /* now compare rest of elements of array to heap */ + for (ia = n; ia < na; ia++) { + /* if new element is greater than smallest, sift it into heap */ + i = 0; + if (ORDER(ia, iord[i])) { + for (ib = (i << 1) + 1; ib < n; ib = (i << 1) + 1) { + if (ib+1 < n) { + if (ORDER(iord[ib], iord[ib+1])) ib++; + } + if (ORDER(ia, iord[ib])) { + iord[i] = iord[ib]; + i = ib; + } else { + break; + } + } + iord[i] = ia; + } + } + + /* unheap iord so largest element is at top */ + for (ia = n - 1; ia > 0; ia--) { + it = iord[ia]; + i = 0; + iord[ia] = iord[i]; + for (ib = (i << 1) + 1; ib < ia; ib = (i << 1) + 1) { + if (ib+1 < ia) { + if (ORDER(iord[ib], iord[ib+1])) ib++; + } + if (ORDER(it, iord[ib])) { + iord[i] = iord[ib]; + i = ib; + } else { + break; + } + } + iord[i] = it; + } +} diff --git a/c_source/vorozobov/readfiles.cpp b/c_source/vorozobov/readfiles.cpp new file mode 100644 index 0000000..d6ac62f --- /dev/null +++ b/c_source/vorozobov/readfiles.cpp @@ -0,0 +1,121 @@ +#include +#include + +#define DL for (d=0;d<3;d++) /* Dimension loop */ +#define BF 1e30 + +/* Positions */ +/* Returns number of particles read */ +int posread(char *posfile, float ***p, float fact) { + + FILE *pos; + int np,dum,d,i; + float xmin,xmax,ymin,ymax,zmin,zmax; + float *ptemp; + + pos = fopen(posfile, "r"); + if (pos == NULL) { + printf("Unable to open position file %s\n\n",posfile); + exit(0); + } + /* Fortran77 4-byte headers and footers */ + /* Delete "dum" statements if you don't need them */ + + /* Read number of particles */ + fread(&dum,1,4,pos); fread(&np,1, sizeof(int),pos); fread(&dum,1,4,pos); + + /* Allocate the arrays */ + (*p) = (float **)malloc(np*sizeof(float *)); + ptemp = (float *)malloc(np*sizeof(float)); + + printf("np = %d\n",np); + + /* Fill the arrays */ + fread(&dum,1,4,pos); + fread(ptemp,np,4,pos); + for (i=0; ixmax) xmax = (*p)[i][0]; + if ((*p)[i][1]ymax) ymax = (*p)[i][1]; + if ((*p)[i][2]zmax) zmax = (*p)[i][2]; + } + printf("np: %d, x: %f,%f; y: %f,%f; z: %f,%f\n",np,xmin,xmax, ymin,ymax, zmin,zmax); fflush(stdout); + + return(np); +} + +/* Velocities */ +/* Returns number of particles read */ +/* Used in voboz, but not zobov, which doesn't use velocities */ +int velread(char *velfile, float ***v, float fact) { + + FILE *vel; + int np,dum,d,i; + float xmin,xmax,ymin,ymax,zmin,zmax; + + vel = fopen(velfile, "r"); + if (vel == NULL) { + printf("Unable to open velocity file %s\n\n",velfile); + exit(0); + } + /* Fortran77 4-byte headers and footers */ + /* Delete "dum" statements if you don't need them */ + + /* Read number of particles */ + fread(&dum,1,4,vel); fread(&np,1, sizeof(int),vel); fread(&dum,1,4,vel); + + /* Allocate the arrays */ + (*v) = (float **)malloc(3*sizeof(float*)); + for (i=0;i<3;i++) (*v)[i] = (float *)malloc(np*sizeof(float)); + + /* Fill the arrays */ + fread(&dum,1,4,vel); fread((*v)[0],np,4,vel); fread(&dum,1,4,vel); + fread(&dum,1,4,vel); fread((*v)[1],np,4,vel); fread(&dum,1,4,vel); + fread(&dum,1,4,vel); fread((*v)[2],np,4,vel); fread(&dum,1,4,vel); + + fclose(vel); + + /* Convert from code units into physical units (km/sec) */ + + for (i=0; i xmax) xmax = (*v)[0][i]; + if ((*v)[1][i] < ymin) ymin = (*v)[1][i]; if ((*v)[1][i] > ymax) ymax = (*v)[1][i]; + if ((*v)[2][i] < zmin) zmin = (*v)[2][i]; if ((*v)[2][i] > zmax) zmax = (*v)[2][i]; + } + printf("vx: %f,%f; vy: %f,%f; vz: %f,%f\n",xmin,xmax, ymin,ymax, zmin,zmax); + + return(np); +} diff --git a/c_source/vorozobov/voz.h b/c_source/vorozobov/voz.h new file mode 100644 index 0000000..f71e1fb --- /dev/null +++ b/c_source/vorozobov/voz.h @@ -0,0 +1,11 @@ +#define MAXVERVER 100000 +#define NGUARD 42 /*Actually, the number of SPACES between guard points + in each dim */ +//#define NGUARD 84 /*Actually, the number of SPACES between guard points + +typedef int pid_t; + +typedef struct Partadj { + pid_t nadj; + pid_t *adj; +} PARTADJ; diff --git a/c_source/vorozobov/voz1b1.cpp b/c_source/vorozobov/voz1b1.cpp new file mode 100644 index 0000000..02574f3 --- /dev/null +++ b/c_source/vorozobov/voz1b1.cpp @@ -0,0 +1,481 @@ +#include "voz.h" +#include +#include +#include +#include + +#include "voro++.hh" +using namespace voro; + +#define DL for (d=0;d<3;d++) +#define BF 1e30 + +#define MAX(a,b) ( ((a) < (b)) ? (a) : (b) ) + +int posread(char *posfile, float ***p, float fact); + +int main(int argc, char *argv[]) { + int exitcode; + int i, j, np; + float **r; + double rtemp[3], *parts; + double deladjs[3*MAXVERVER], points[3*MAXVERVER]; + FILE *pos, *out, *testFile; + char *posfile, outfile[200], *suffix, *outDir; + PARTADJ *adjs; + float *vols; + float predict, xmin,xmax,ymin,ymax,zmin,zmax; + int *orig; + + int isitinbuf; + char isitinmain, d; + int numdiv; + int nvp, nvpall, nvpbuf; + float width, width2, totwidth, totwidth2, bf, s, g; + float border, boxsize; + float c[3]; + int b[3]; + double totalvol; + int isObs = 0; + + printf("Routine: voz1b1\n"); + + if (argc != 11) { + printf("Wrong number of arguments.\n"); + printf("arg1: position file\n"); + printf("arg2: border size\n"); + printf("arg3: boxsize\n"); + printf("arg4: suffix\n"); + printf("arg5: number of divisions\n"); + printf("arg6-8: b[0-2]\n\n"); + printf("arg9: output directory\n"); + printf("arg10: observation flag\n"); + exit(0); + } + posfile = argv[1]; + if (sscanf(argv[2],"%f",&border) != 1) { + printf("That's no border size; try again.\n"); + exit(0); + } + if (sscanf(argv[3],"%f",&boxsize) != 1) { + printf("That's no boxsize; try again.\n"); + exit(0); + } + suffix = argv[4]; + if (sscanf(argv[5],"%d",&numdiv) != 1) { + printf("%s is no number of divisions; try again.\n",argv[5]); + exit(0); + } + if (numdiv == 1) { + printf("Only using one division; should only use for an isolated segment.\n"); + } + if (numdiv < 1) { + printf("Cannot have a number of divisions less than 1. Resetting to 1.\n"); + numdiv = 1; + } + if (sscanf(argv[6],"%d",&b[0]) != 1) { + printf("That's no b index; try again.\n"); + exit(0); + } + if (sscanf(argv[7],"%d",&b[1]) != 1) { + printf("That's no b index; try again.\n"); + exit(0); + } + if (sscanf(argv[8],"%d",&b[2]) != 1) { + printf("That's no b index; try again.\n"); + exit(0); + } + outDir = argv[9]; + if (sscanf(argv[10],"%d",&isObs) != 1) { + printf("That's no observation mode; try again.\n"); + exit(0); + } + + printf("Working with subdomain b=[%d,%d,%d]\n", b[0], b[1], b[2]); + + /* Boxsize should be the range in r, yielding a range 0-1 */ + printf("Reading particle file...\n"); + np = posread(posfile,&r,1./boxsize); + printf("%d particles\n",np);fflush(stdout); + + xmin = BF; xmax = -BF; ymin = BF; ymax = -BF; zmin = BF; zmax = -BF; + for (i=0; ixmax) xmax = r[i][0]; + if (r[i][1]ymax) ymax = r[i][1]; + if (r[i][2]zmax) zmax = r[i][2]; + } + printf("Full box particle stats: np: %d, x: %f,%f; y: %f,%f; z: %f,%f\n",np, + xmin,xmax, ymin,ymax, zmin,zmax); fflush(stdout); + + width = 1./(float)numdiv; + width2 = 0.5*width; + if (border > 0.) bf = border; + else bf = 0.1; + /* In units of 0-1, the thickness of each subregion's buffer*/ + totwidth = width+2.*bf; + totwidth2 = width2 + bf; + + s = width/(float)NGUARD; + if ((bf*bf - 2.*s*s) < 0.) { + printf("bf = %f, s = %f.\n",bf,s); + printf("Not enough guard points for given border.\nIncrease guards to >= %f\n.", + sqrt(2.)*width/bf); + exit(0); + } + g = (bf / 2.)*(1. + sqrt(1 - 2.*s*s/(bf*bf))); + printf("Guard particle information: s = %f, bf = %f, g = %f.\n",s,bf,g); + + fflush(stdout); + + adjs = (PARTADJ *)malloc(np*sizeof(PARTADJ)); + if (adjs == NULL) { + printf("Unable to allocate adjs\n"); + exit(1); + } + + DL c[d] = ((float)b[d]+0.5)*width; + printf("Counting particles in subdomain...\n"); + /* Assign temporary array*/ + nvpbuf = 0; /* Number of particles to tesselate, includng buffer and guards */ + nvp = 0; /* Without the buffer */ + for (i=0; i 0.5) rtemp[d] --; + if (rtemp[d] < -0.5) rtemp[d] ++; + isitinbuf = isitinbuf && (fabs(rtemp[d]) < totwidth2); + isitinmain = isitinmain && (fabs(rtemp[d]) <= width2); + } + + if (isitinbuf) nvpbuf++; + if (isitinmain) nvp++; + } + + nvpbuf += 6*(NGUARD+1)*(NGUARD+1); /* number of guard points */ + + parts = (double *)malloc(3*nvpbuf*sizeof(double)); + orig = (pid_t *)malloc(nvpbuf*sizeof(pid_t)); + + if (parts == NULL) { + printf("Unable to allocate parts\n"); + fflush(stdout); + exit(1); + } + if (orig == NULL) { + printf("Unable to allocate orig\n"); + fflush(stdout); + exit(1); + } + + printf("Placing particles in subdomain...\n"); + nvp = 0; nvpall = 0; /* nvp = number of particles without buffer */ + xmin = BF; xmax = -BF; ymin = BF; ymax = -BF; zmin = BF; zmax = -BF; + for (i=0; i 0.5) rtemp[d] --; + if (rtemp[d] < -0.5) rtemp[d] ++; + isitinmain = isitinmain && (fabs(rtemp[d]) <= width2); + } + if (isitinmain) { + parts[3*nvp] = rtemp[0]; + parts[3*nvp+1] = rtemp[1]; + parts[3*nvp+2] = rtemp[2]; + orig[nvp] = i; + nvp++; + if (rtemp[0] < xmin) xmin = rtemp[0]; + if (rtemp[0] > xmax) xmax = rtemp[0]; + if (rtemp[1] < ymin) ymin = rtemp[1]; + if (rtemp[1] > ymax) ymax = rtemp[1]; + if (rtemp[2] < zmin) zmin = rtemp[2]; + if (rtemp[2] > zmax) zmax = rtemp[2]; + } + } + printf("Subdomain: b=(%d,%d,%d), c=(%f,%f,%f)\n", + b[0],b[1],b[2],c[0],c[1],c[2]); + printf(" nvp=%d, nvpbuf=%d\n", nvp,nvpbuf); + printf(" Particle ranges: x: %f,%f; y: %f,%f; z:%f,%f\n", + xmin,xmax,ymin,ymax,zmin,zmax); + + printf("Adding particles to buffer regions...\n"); + nvpbuf = nvp; + for (i=0; i 0.5) rtemp[d] --; + if (rtemp[d] < -0.5) rtemp[d] ++; + isitinbuf = isitinbuf && (fabs(rtemp[d]) xmax) xmax = rtemp[0]; + if (rtemp[1] < ymin) ymin = rtemp[1]; + if (rtemp[1] > ymax) ymax = rtemp[1]; + if (rtemp[2] < zmin) zmin = rtemp[2]; + if (rtemp[2] > zmax) zmax = rtemp[2]; + } + } + nvpall = nvpbuf; + predict = pow(width+2.*bf,3)*(float)np; + printf("Total particles including buffer = %d (predicted ~%f)\n", + nvpbuf, predict); + printf(" Particle ranges including buffer: x: %f,%f; y: %f,%f; z:%f,%f\n", + xmin,xmax,ymin,ymax,zmin,zmax); + + for (i=0;i nvp) itype = 1; // it's a buffer particle + if (i > nvpbuf) itype = 2; // it's a guard particle + fprintf(testFile, "%f %f %f %d %f %f %f\n", c[0], c[1], c[2], + itype, + parts[3*i], + parts[3*i+1], + parts[3*i+2]); + } + fclose(testFile); + + 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]; + } + + printf("Added guard points to total %d points (should be %d)\n",nvpall, + nvpbuf + 6*(NGUARD+1)*(NGUARD+1)); + printf("New particle ranges including guards: x: %f,%f; y: %f,%f; z:%f,%f\n", + xmin,xmax,ymin,ymax,zmin,zmax); + + /* Do tesselation*/ + printf("File read. Tessellating ...\n"); fflush(stdout); + + // initialize a voro++ container with non-periodic boundaries + // per voro++ documentaiton, max efficiency is achieved with ~8 particles per subdomain block + int nblocks = floor(pow(nvpall, 1/3.)/8); + container con(xmin,xmax,ymin,ymax,zmin,zmax, nblocks,nblocks,nblocks, false,false,false, 8); + + particle_order po; // this will keep track of only original particles + + // initialize storage arrays + vols = (float *)malloc(nvp*sizeof(float)); + + // place our particles in the container + for (i=0; i adjsv; + cell.neighbors(adjsv); + + // copy adjacencies to old-school data format + adjs[p].nadj = adjsv.size(); + adjs[p].adj = (int *)malloc(adjsv.size() * sizeof(int)); + + for (int iAdj; iAdj < adjsv.size(); ++iAdj) { + adjs[p].adj[iAdj] = adjsv[iAdj]; + } + + // grab volume and store + vols[p] = cell.volume()*(float)np; + p++; + } while(vl.inc()); + + //if (i % 1000 == 0) + // printf("%d: %d, %f\n",i,adjs[i].nadj,vols[i]); + +// PMS - reset number of adjancies to not include links to border guards +/* + for (i=0; i 0) + fwrite(adjs[i].adj,adjs[i].nadj,sizeof(pid_t),out); + else printf("0"); + } + fclose(out); + + return(0); +} diff --git a/c_source/vorozobov/vozinit.cpp b/c_source/vorozobov/vozinit.cpp new file mode 100644 index 0000000..5a3fe58 --- /dev/null +++ b/c_source/vorozobov/vozinit.cpp @@ -0,0 +1,203 @@ +#include "voz.h" +#include +#include +#include + +#define NUMCPU 2 +#define DL for (d=0;d<3;d++) +#define BF 1e30 + +#define QHULL_MAX_PARTICLES ((1L<<24)-1) + +int posread(char *posfile, float ***p, float fact); + +int main(int argc, char *argv[]) { + int i, np; + float **rfloat, rtemp[3]; + FILE *pos, *scr; + char *posfile, scrfile[200], systemstr[90], *suffix, *outDir, *vobozPath; + float xmin,xmax,ymin,ymax,zmin,zmax; + + int isitinbuf; + char isitinmain, d; + int numdiv; + int p; + int nvp, nvpall, nvpbuf, nvpmin, nvpmax, nvpbufmin, nvpbufmax; /* yes, the insurance */ + float width, width2, totwidth, totwidth2, bf, s, g; + float widthX, widthY, widthZ; + float border, boxsize, boxsizeX, boxsizeY, boxsizeZ; + float c[3]; + int numGuards; + int b[3]; + int numThreads; + int isObs = 0; + + printf("Routine: vozinit (voro++ implementation)\n"); + + if (argc != 10) { + printf("Wrong number of arguments.\n"); + printf("arg1: position file\n"); + printf("arg2: buffer size (default 0.1)\n"); + printf("arg3: box size\n"); + printf("arg4: number of divisions (default 2)\n"); + printf("arg5: suffix describing this run\n"); + printf("arg6: number of parallel threads\n"); + printf("arg7: location of voboz executables\n"); + printf("arg8: output directory\n"); + printf("arg9: observation flag (set to 1 for observation mode)\n"); + exit(0); + } + posfile = argv[1]; + suffix = argv[2]; + if (sscanf(suffix,"%f",&border) != 1) { + printf("That's no border size; try again.\n"); + exit(0); + } + suffix = argv[3]; + if (sscanf(suffix,"%f",&boxsize) != 1) { + printf("That's no boxsize; try again.\n"); + exit(0); + } + suffix = argv[4]; + if (sscanf(suffix,"%d",&numdiv) != 1) { + printf("That's no number of divisions; try again.\n"); + exit(0); + } + if (numdiv < 1) { + printf("Cannot have a number of divisions less than 1. Resetting to 1:\n"); + numdiv = 1; + } + + suffix = argv[5]; + vobozPath = argv[6]; + if (sscanf(vobozPath,"%d",&numThreads) != 1) { + printf("That's no number of threads; try again.\n"); + exit(0); + } + vobozPath = argv[7]; + outDir = argv[8]; + if (sscanf(argv[9],"%d",&isObs) != 1) { + printf("That's no observation mode; try again.\n"); + exit(0); + } + + /* Read the position file */ + printf("Reading particle file and calculating box sizes...\n"); + np = posread(posfile,&rfloat,1./boxsize); + /* Boxsize should be the range in r, yielding a range 0-1 */ + + width = boxsize/(float)numdiv; + //width = 1./(float)numdiv; + width2 = 0.5*width; + if (border > 0.) bf = border; + else bf = 0.1; + + /* In units of 0-1, the thickness of each subregion's buffer*/ + totwidth = width+2.*bf; + totwidth2 = width2 + bf; + + s = width/(float)NGUARD; + if ((bf*bf - 2.*s*s) < 0.) { + printf("Not enough guard points for given border.\nIncrease guards to >= %f\n.", + totwidth/sqrt(0.5*bf*bf)); + printf("bf = %f\n",bf); + exit(0); + } + g = (bf / 2.)*(1. + sqrt(1 - 2.*s*s/(bf*bf))); + printf("s = %f, bf = %f, g = %f.\n",s,bf,g); + + nvpmax = 0; nvpbufmax = 0; nvpmin = np; nvpbufmin = np; + + // count particles in each sub-domain + printf("Counting particles in each subdomain...\n"); + for (b[0] = 0; b[0] < numdiv; b[0]++) { + c[0] = ((float)b[0]+0.5)*width; + for (b[1] = 0; b[1] < numdiv; b[1]++) { + c[1] = ((float)b[1]+0.5)*width; + for (b[2] = 0; b[2] < numdiv; b[2]++) { + c[2] = ((float)b[2]+0.5)*width; + + nvp = 0; /* Number of particles excluding buffer */ + nvpbuf = 0; /* Number of particles to tesselate, including buffer */ + xmin = BF; xmax = -BF; ymin = BF; ymax = -BF; zmin = BF; zmax = -BF; + + for (i=0; i 0.5) rtemp[d] --; + if (rtemp[d] < -0.5) rtemp[d] ++; + //} + isitinbuf = isitinbuf && (fabs(rtemp[d]) < totwidth2); + isitinmain = isitinmain && (fabs(rtemp[d]) <= width2); + } + if (isitinbuf) { nvpbuf++; } + if (isitinmain) { + nvp++; + if (rtemp[0] < xmin) xmin = rtemp[0]; + if (rtemp[0] > xmax) xmax = rtemp[0]; + if (rtemp[1] < ymin) ymin = rtemp[1]; + if (rtemp[1] > ymax) ymax = rtemp[1]; + if (rtemp[2] < zmin) zmin = rtemp[2]; + if (rtemp[2] > zmax) zmax = rtemp[2]; + } + } + + if (nvp > nvpmax) nvpmax = nvp; + if (nvpbuf > nvpbufmax) nvpbufmax = nvpbuf; + if (nvp < nvpmin) nvpmin = nvp; + if (nvpbuf < nvpbufmin) nvpbufmin = nvpbuf; + + printf("Subdomain: b=(%d,%d,%d), c=(%f,%f,%f)\n", + b[0],b[1],b[2],c[0],c[1],c[2]); + printf(" nvp=%d, nvpbuf=%d\n", nvp,nvpbuf); + printf(" Particle ranges: x: %f,%f; y: %f,%f; z:%f,%f\n", + xmin,xmax,ymin,ymax,zmin,zmax); + } + } + } + printf("Nvp range: %d,%d\n",nvpmin,nvpmax); + printf("Nvpbuf range: %d,%d\n",nvpbufmin,nvpbufmax); + + numGuards = 6*(NGUARD+1)*(NGUARD+1); + printf("Num guards = %d\n", numGuards); + printf("Total max particles: %d\n" , nvpbufmax+numGuards); + if (nvpbufmax+numGuards >= QHULL_MAX_PARTICLES) { + printf("Too many particles to tesselate per division (Qhull will overflow). Increase divisions or reduce buffer size.\n"); + fflush(stdout); + exit(1); + } + + /* Output script file */ + sprintf(scrfile,"scr%s",suffix); + printf("Writing script file to %s.\n",scrfile);fflush(stdout); + scr = fopen(scrfile,"w"); + if (scr == NULL) { + printf("Problem opening script file.\n"); + fflush(stdout); + exit(1); + } + fprintf(scr,"#!/bin/bash -f\n"); + p = 0; + for (b[0]=0;b[0] +#include +#include +#include +#include +#include "voz.h" + +int compare_int(const void *a, const void *b) +{ + const int *ia = (const int *)a; + const int *ib = (const int *)b; + if ((*ia) < (*ib)) + return -1; + else if ((*ia) > (*ib)) + return 1; + else + return 0; +} + +int main(int argc, char *argv[]) { + + FILE *part, *adj, *vol; + char partfile[200], *suffix, adjfile[200], volfile[200], *outDir; + float *vols, volstemp; + int *cnt_adj; + + PARTADJ *adjs; + + int numdiv,np,np2,na; + + pid_t i,j,k,p,nout; + pid_t nvp,npnotdone,nvpmax, nvpsum, *orig; + double avgnadj, avgvol; + + int numRemoved = 0; + int isObs = 0; + + if (argc != 5) { + printf("Voztie: Wrong number of arguments.\n"); + printf("arg1: number of divisions (default 2)\n"); + printf("arg2: suffix describing this run\n"); + printf("arg3: file directory\n"); + printf("arg5: observation flag\n"); + exit(0); + } + if (sscanf(argv[1],"%d",&numdiv) != 1) { + printf("That's no number of divisions; try again.\n"); + exit(0); + } + if (numdiv < 1) { + printf("Cannot have a number of divisions less than 1. Resetting to 1:\n"); + numdiv = 1; + } + + suffix = argv[2]; + outDir = argv[3]; + + if (sscanf(argv[4],"%d",&isObs) != 1) { + printf("That's no observation flag; try again.\n"); + exit(0); + } + + np = -1; nvpmax = -1; nvpsum = 0; + + for (i = 0; i < numdiv; i++) { + for (j = 0; j < numdiv; j++) { + for (k = 0; k < numdiv; k++) { + sprintf(partfile,"%s/part.%s.%02d.%02d.%02d",outDir,suffix,i,j,k); + part = fopen(partfile,"r"); + if (part == NULL) { + printf("Unable to open file %s.\n\n",partfile); + exit(0); + } + fread(&np2,1,sizeof(int),part); + fread(&nvp,1,sizeof(int),part); + nvpsum += nvp; + if (np == -1) + np = np2; + else + if (np2 != np) { + printf("Incompatible total particle numbers: %d,%d\n\n",np,np2); + exit(0); + } + if (nvp > nvpmax) nvpmax = nvp; + fclose(part); + } + } + } + + printf("We have %d particles to tie together.\n",np); fflush(stdout); + printf("The maximum number of particles in a file is %d.\n",nvpmax); + + cnt_adj = (int *)malloc(np*sizeof(int)); + if (cnt_adj == NULL) + printf("Could not allocate cnt_adj.\n"); + + adjs = (PARTADJ *)malloc(np*sizeof(PARTADJ)); + if (adjs == NULL) printf("Couldn't allocate adjs.\n"); + vols = (float *)malloc(np*sizeof(float)); + if (vols == NULL) printf("Couldn't allocate vols.\n"); + orig = (pid_t *)malloc(nvpmax*sizeof(pid_t)); + if (orig == NULL) printf("Couldn't allocate orig.\n"); + if ((cnt_adj==NULL) || (vols == NULL) || (orig == NULL) || (adjs == NULL)) { + printf("Not enough memory to allocate. Exiting.\n"); + exit(0); + } + for (p=0;p -1.) + if (fabs(vols[orig[p]]-volstemp)/volstemp > 1.5e-3) { + printf("Inconsistent volumes for p. %d: (%10g,%10g)!\n", + orig[p],vols[orig[p]],volstemp); + //exit(0); + } + vols[orig[p]] = volstemp; + } + + for (p=0;p 0) { + pid_t *previous_adj = adjs[pid].adj; + assert(adjs[pid].nadj == 0 || adjs[pid].nadj == na); + adjs[pid].nadj = na; + adjs[pid].adj = (pid_t *)malloc(na*sizeof(pid_t)); + if (adjs[pid].adj == 0) { + printf("Couldn't allocate adjs[orig[%d]].adj.\n",p); + exit(0); + } + fread(adjs[pid].adj,na,sizeof(pid_t),part); + if (previous_adj != 0) + { + qsort(previous_adj, na, sizeof(pid_t), compare_int); + qsort(adjs[pid].adj, na, sizeof(pid_t), compare_int); + if (memcmp(previous_adj, adjs[pid].adj, sizeof(pid_t)*na) != 0) + { + printf("Inconsistent adjacencies between two divisions.\n"); + abort(); + } + } + } else { + printf("0"); fflush(stdout); + } + } + fclose(part); + printf("%d ",k); + } + } + } + printf("\n"); + +/* + // PMS : remove mock galaxies and anything adjacent to a mock galaxy + printf("\nRemoving mock galaxies...\n"); + + // unlink particles adjacent to mock galaxies + for (i = 0; i < np; i++) { + for (j = 0; j < adjs[i].nadj; j++) { + if (adjs[i].adj[j] > np) { +//printf("KILLING %d\n", i); + vols[i] = 1.e-29; + adjs[i].nadj = 0; + numRemoved++; + break; + } + } + } + + // update all other adjacencies + for (i = 0; i < np; i++) { + + int numAdjSaved = 0; + for (j = 0; j < adjs[i].nadj; j++) { + + if ( adjs[adjs[i].adj[j]].nadj != 0) { + adjs[i].adj[numAdjSaved] = adjs[i].adj[j]; + numAdjSaved++; + } + + } + adjs[i].nadj = numAdjSaved; + } + + // END PMS +*/ + + npnotdone = 0; avgnadj = 0.; avgvol = 0.; + for (p=0;p 0) + printf("%d particles not done!\n", npnotdone); + printf("%d particles done more than once.\n",nvpsum-np); + avgnadj /= (double)np; + avgvol /= (double)np; + printf("Average # adjacencies = %lf (%f for Poisson)\n",avgnadj, + 48.*3.141593*3.141593/35.+2.); + printf("Average volume = %lf\n",avgvol); + + /* Now the output! */ + + sprintf(adjfile,"%s/adj%s.dat",outDir,suffix); + sprintf(volfile,"%s/vol%s.dat",outDir,suffix); + + printf("Outputting to%s, %s\n\n", adjfile, volfile); + + adj = fopen(adjfile,"w"); + if (adj == NULL) { + printf("Unable to open %s\n",adjfile); + exit(0); + } + fwrite(&np,1, sizeof(int),adj); + /* Adjacencies: first the numbers of adjacencies, + and the number we're actually going to write per particle */ + + // Recount the number of adjacencies after merge + for(i=0;i i) + nout++; + fwrite(&nout,1,sizeof(int),adj); + for (j=0;j i) + fwrite(&id,1,sizeof(pid_t),adj); + } + } + + fclose(adj); + + /* Volumes */ + vol = fopen(volfile,"w"); + if (vol == NULL) { + printf("Unable to open %s\n",volfile); + exit(0); + } + fwrite(&np,1, sizeof(int),vol); + fwrite(vols,sizeof(float),np,vol); + + fclose(vol); + + free(vols); + free(cnt_adj); + free(adjs); + free(orig); + + return(0); +} diff --git a/c_source/vorozobov/vozutil.cpp b/c_source/vorozobov/vozutil.cpp new file mode 100644 index 0000000..7eb9bae --- /dev/null +++ b/c_source/vorozobov/vozutil.cpp @@ -0,0 +1,250 @@ +#include "voz.h" +#include +#include +#include + +#include "voro++.hh" +using namespace voro; + +#define FOREACHvertex2_(vertices) FOREACHsetelement_(vertexT, vertices2,vertex2) + +int compar(const void * n1, const void * n2) { + int i1,i2; + + i1 = *(int *)n1; + i2 = *(int *)n2; + return 2*(i1 > i2) - 1 + (i1 == i2); +} + +/* Finds Delaunay adjacencies of a set of points */ +int delaunadj (double *points, int nvp, int nvpbuf, int nvpall, PARTADJ **adjs) { + + container con(0.0,1.0,0.0,1.0,0.0,1.0,10,10,10,false,false,false,8); + + int dim= 3; /* dimension of points */ + bool ismalloc= False; /* True if qhull should free points in qh_freeqhull() or reallocation */ + char flags[250]; /* option flags for qhull, see qh_opt.htm */ + FILE *outfile= stdout; /* output from qh_produce_output() + use NULL to skip qh_produce_output() */ + FILE *errfile= stderr; /* error messages from qhull code */ + int exitcode; /* 0 if no error from qhull */ + int curlong, totlong; /* memory remaining after qh_memfreeshort */ + int i, ver, count; + int numfacets, numsimplicial, numridges, totneighbors, numneighbors, + numcoplanars, numtricoplanars; + setT *vertices, *vertices2, *vertex_points, *coplanar_points; + vertexT *vertex, **vertexp; + vertexT *vertex2, **vertex2p; + int vertex_i, vertex_n; + facetT *facet, **facetp, *neighbor, **neighborp; + pointT *point, **pointp; + int numdiv; + + PARTADJ adjst; + + int errorReported = 0; + + adjst.adj = (int *)malloc(MAXVERVER*sizeof(int)); + if (adjst.adj == NULL) { + printf("Unable to allocate adjst.adj\n"); + exit(0); + } + + /* Delaunay triangulation*/ + sprintf (flags, "qhull s d Qt"); + exitcode= qh_new_qhull (dim, nvpall, points, ismalloc, + flags, outfile, errfile); + + + if (!exitcode) { /* if no error */ + /* 'qh facet_list' contains the convex hull */ + + /* From qh_printvneighbors */ + qh_countfacets(qh facet_list, NULL, 0, &numfacets, &numsimplicial, + &totneighbors, &numridges, &numcoplanars, &numtricoplanars); + qh_vertexneighbors(); + vertices= qh_facetvertices (qh facet_list, NULL, 0); + vertex_points= qh_settemp (nvpall); + coplanar_points= qh_settemp (nvpall); + qh_setzero (vertex_points, 0, nvpall); + qh_setzero (coplanar_points, 0, nvpall); + FOREACHvertex_(vertices) + qh_point_add (vertex_points, vertex->point, vertex); + FORALLfacet_(qh facet_list) { + FOREACHpoint_(facet->coplanarset) + qh_point_add (coplanar_points, point, facet); + } + ver = 0; + FOREACHvertex_i_(vertex_points) { + (*adjs)[ver].nadj = 0; + if (vertex) { + /* Count the neighboring vertices, check that all are real + neighbors */ + adjst.nadj = 0; + FOREACHneighbor_(vertex) { + if ((*adjs)[ver].nadj > -1) { + if (neighbor->visitid) { + vertices2 = neighbor->vertices; + FOREACHvertex2_(vertices2) { + if (ver != qh_pointid(vertex2->point)) { + adjst.adj[adjst.nadj] = (int)qh_pointid(vertex2->point); + adjst.nadj ++; + if (adjst.nadj > MAXVERVER) { + printf("Increase MAXVERVER to at least %d!\n",adjst.nadj); + exit(0); + } + } + } + } else { + printf(" %d",ver); + (*adjs)[ver].nadj = -1; /* There are unreal vertices here */ + } + } + } + } else (*adjs)[ver].nadj = -2; + + /* Enumerate the unique adjacencies*/ + if (adjst.nadj >= 4) { + qsort((void *)adjst.adj, adjst.nadj, sizeof(int), &compar); + count = 1; + + for (i=1; i= nvpbuf && !errorReported) { + errorReported = 1; + printf("Guard point encountered. Increase border and/or nguard.\n"); + printf("P:(%f,%f,%f), G: (%f,%f,%f)\n",points[3*ver],points[3*ver+1],points[3*ver+2], + points[3*adjst.adj[i]],points[3*adjst.adj[i]+1],points[3*adjst.adj[i]+2]); + } + count++; + } + (*adjs)[ver].adj = (int *)malloc(count*sizeof(int)); + if ((*adjs)[ver].adj == NULL) { + printf("Unable to allocate (*adjs)[ver].adj\n"); + exit(0); + } + (*adjs)[ver].adj[0] = adjst.adj[0]; + count = 1; + for (i=1; i %d\n",adjst.nadj,ver,ver); + exit(0); + } + ver++; + if (ver == nvp) break; + } + qh_settempfree (&coplanar_points); + qh_settempfree (&vertex_points); + qh_settempfree (&vertices); + } + qh_freeqhull(!qh_ALL); /* free long memory */ + qh_memfreeshort (&curlong, &totlong); /* free short memory and memory allocator */ + if (curlong || totlong) + fprintf (errfile, "qhull internal warning (delaunadj): did not free %d bytes of long memory (%d pieces)\n", totlong, curlong); + free(adjst.adj); + return exitcode; +} + +/* Calculates the Voronoi volume from a set of Delaunay adjacencies */ +int vorvol (double *deladjs, double *points, pointT *intpoints, int numpoints, float *vol) { + int dim= 3; /* dimension of points */ + bool ismalloc= False; /* True if qhull should free points in qh_freeqhull() or reallocation */ + char flags[250]; /* option flags for qhull, see qh_opt.htm */ + FILE *outfile= NULL; /* output from qh_produce_output() + use NULL to skip qh_produce_output() */ + FILE *errfile= stderr; /* error messages from qhull code */ + int exitcode; /* 0 if no error from qhull */ + facetT *facet; /* set by FORALLfacets */ + int curlong, totlong; /* memory remaining after qh_memfreeshort */ + + double *point, *normp, *coordp, *feasiblep, *deladj; + int i, j, k; + bool zerodiv; + float runsum; + char region; + /*double *points; + pointT *intpoints;*/ + + /* make point array from adjacency coordinates (add offset)*/ + /*points = (double *)malloc(4*numpoints*sizeof(double)); + if (points == NULL) { + printf("Unable to allocate points\n"); + exit(0); + }*/ + for (i=0; inormal; + feasiblep= qh feasible_point; + if (facet->offset < -qh MINdenom) { + for (k= qh hull_dim; k--; ) + *(coordp++)= (*(normp++) / - facet->offset) + *(feasiblep++); + }else { + for (k= qh hull_dim; k--; ) { + *(coordp++)= qh_divzero (*(normp++), facet->offset, qh MINdenom_1, + &zerodiv) + *(feasiblep++); + if (zerodiv) { + qh_memfree (point, qh normal_size); + printf("LABELprintinfinite\n"); + exit(0); + } + } + } + } + } + qh_freeqhull (!qh_ALL); + qh_memfreeshort (&curlong, &totlong); + + /* Now we calculate the volume */ + sprintf (flags, "qhull FA"); + exitcode= qh_new_qhull (dim, numpoints, intpoints, ismalloc, + flags, outfile, errfile); + + qh_getarea(qh facet_list); + *vol = qh totvol; + + qh_freeqhull (!qh_ALL); + qh_memfreeshort (&curlong, &totlong); + if (curlong || totlong) + fprintf (errfile, "qhull internal warning (vorvol): did not free %d bytes of long memory (%d pieces)\n", totlong, curlong); + /*free(points); free(intpoints);*/ + + return exitcode; +} diff --git a/external/external_build.cmake b/external/external_build.cmake index 0abb8f7..a0f8b4d 100644 --- a/external/external_build.cmake +++ b/external/external_build.cmake @@ -398,4 +398,5 @@ include_directories(${CMAKE_BINARY_DIR}/src ${Boost_INCLUDE_DIRS} ${QHULL_INCLUDE_PATH} ${LIBSDF_INCLUDE_PATH}) + message(STATUS "Boost_INTERNAL_INSTALL=${Boost_INTERNAL_INSTALL}")