mirror of
https://bitbucket.org/cosmicvoids/vide_public.git
synced 2025-07-05 07:41:11 +00:00
Initial import
This commit is contained in:
Guilhem Lavaux
commit
e57d9a731a
17 changed files with 3271 additions and 0 deletions
7
.gitignore
vendored
Normal file
7
.gitignore
vendored
Normal file
|
|
@ -0,0 +1,7 @@
|
||||||
|
*~
|
||||||
|
*.o
|
||||||
|
*.prog
|
||||||
|
voz1b1
|
||||||
|
jozov
|
||||||
|
vozinit
|
||||||
|
voztie
|
||||||
41
Makefile
Normal file
41
Makefile
Normal file
|
|
@ -0,0 +1,41 @@
|
||||||
|
# Compiler choice:
|
||||||
|
|
||||||
|
# Gcc
|
||||||
|
CC = icc
|
||||||
|
CFLAGS = -O3 -ggdb
|
||||||
|
|
||||||
|
MLIBS = -lm
|
||||||
|
|
||||||
|
#QLIB = -L../qhull-2003.1/src/
|
||||||
|
#QINC = -I../qhull-2003.1/src/
|
||||||
|
QLIB = -L../qhull2002.1/src/
|
||||||
|
QINC = -I../qhull2002.1/src/
|
||||||
|
|
||||||
|
###############
|
||||||
|
|
||||||
|
all: vozinit voz1b1 voztie jozov
|
||||||
|
|
||||||
|
jozov: jozov.o findrtop.o Makefile
|
||||||
|
$(CC) $(CFLAGS) -o jozov jozov.o findrtop.o $(MLIBS)
|
||||||
|
jozov.o: jozov.c Makefile
|
||||||
|
$(CC) $(CFLAGS) -c jozov.c
|
||||||
|
findrtop.o: findrtop.c Makefile
|
||||||
|
$(CC) $(CFLAGS) -c findrtop.c
|
||||||
|
|
||||||
|
voz1b1: voz1b1.o readfiles.o vozutil.o voz.h Makefile
|
||||||
|
$(CC) -o voz1b1 $(CFLAGS) voz1b1.o readfiles.o vozutil.o -L. $(QLIB) -lqhull $(MLIBS)
|
||||||
|
voz1b1.o: voz1b1.c Makefile
|
||||||
|
$(CC) $(CFLAGS) $(QINC) -c voz1b1.c
|
||||||
|
vozutil.o: vozutil.c Makefile
|
||||||
|
$(CC) $(CFLAGS) $(QINC) -c vozutil.c
|
||||||
|
|
||||||
|
vozinit: vozinit.o readfiles.o voz.h Makefile
|
||||||
|
$(CC) -o vozinit $(CFLAGS) vozinit.o readfiles.o -L. $(MLIBS)
|
||||||
|
vozinit.o: vozinit.c Makefile
|
||||||
|
$(CC) $(CFLAGS) -c vozinit.c $(QINC)
|
||||||
|
|
||||||
|
voztie: voztie.o readfiles.o Makefile
|
||||||
|
$(CC) -o voztie $(CFLAGS) voztie.o readfiles.o
|
||||||
|
voztie.o: voztie.c Makefile
|
||||||
|
$(CC) $(CFLAGS) -c voztie.c
|
||||||
|
|
||||||
81
findrtop.c
Normal file
81
findrtop.c
Normal file
|
|
@ -0,0 +1,81 @@
|
||||||
|
#include <math.h>
|
||||||
|
#include <stdlib.h>
|
||||||
|
|
||||||
|
/*------------------------------------------------------------------------------
|
||||||
|
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;
|
||||||
|
}
|
||||||
|
}
|
||||||
550
jozov.c
Normal file
550
jozov.c
Normal file
|
|
@ -0,0 +1,550 @@
|
||||||
|
#include <assert.h>
|
||||||
|
/* jovoz.c by Mark Neyrinck */
|
||||||
|
#include <stdio.h>
|
||||||
|
#include <stdlib.h>
|
||||||
|
#include <math.h>
|
||||||
|
|
||||||
|
#define BIGFLT 1e30 /* Biggest possible floating-point number */
|
||||||
|
#define NLINKS 100 /* Number of possible links with the same rho_sl */
|
||||||
|
#define FF fflush(stdout)
|
||||||
|
|
||||||
|
typedef struct Particle {
|
||||||
|
float dens;
|
||||||
|
int nadj;
|
||||||
|
int ncnt;
|
||||||
|
int *adj;
|
||||||
|
} PARTICLE;
|
||||||
|
|
||||||
|
typedef struct Zone {
|
||||||
|
int core; /* Identity of peak particle */
|
||||||
|
int np; /* Number of particles in zone */
|
||||||
|
int npjoin; /* Number of particles in the joined void */
|
||||||
|
int nadj; /* Number of adjacent zones */
|
||||||
|
int nhl; /* Number of zones in final joined void */
|
||||||
|
float leak; /* Volume of last leak zone*/
|
||||||
|
int *adj; /* Each adjacent zone, with ... */
|
||||||
|
float *slv; /* Smallest Linking Volume */
|
||||||
|
float denscontrast; /* density contrast */
|
||||||
|
double vol; /* Total volume of all particles in the zone */
|
||||||
|
double voljoin; /* Total volume of all particles in the joined void */
|
||||||
|
} ZONE;
|
||||||
|
|
||||||
|
typedef struct ZoneT {
|
||||||
|
int nadj; /* Number of zones on border */
|
||||||
|
int *adj; /* Each adjacent zone, with ... */
|
||||||
|
float *slv; /* Smallest Linking Volume */
|
||||||
|
} ZONET;
|
||||||
|
|
||||||
|
void findrtop(double *a, int na, int *iord, int nb);
|
||||||
|
|
||||||
|
int main(int argc,char **argv) {
|
||||||
|
|
||||||
|
FILE *adj, *vol, *zon, *zon2, *txt;
|
||||||
|
PARTICLE *p;
|
||||||
|
ZONE *z;
|
||||||
|
ZONET *zt;
|
||||||
|
char *adjfile, *volfile, *zonfile, *zonfile2, *txtfile;
|
||||||
|
int i, j,k,l, h, h2,hl,n,np, np2, nzones, nhl, nhlcount, nhl2;
|
||||||
|
int *jumper, *jumped, *numinh;
|
||||||
|
int *zonenum, *zonelist, *zonelist2;
|
||||||
|
int link[NLINKS], link2, nl;
|
||||||
|
float lowvol, voltol, prob;
|
||||||
|
|
||||||
|
int q,q2;
|
||||||
|
|
||||||
|
int za, nin;
|
||||||
|
int testpart;
|
||||||
|
char already, interior, *inyet, *inyet2, added, beaten;
|
||||||
|
int *nm, **m;
|
||||||
|
|
||||||
|
float maxvol, minvol;
|
||||||
|
double *sorter, e1,maxdenscontrast;
|
||||||
|
int *iord;
|
||||||
|
|
||||||
|
e1 = exp(1.)-1.;
|
||||||
|
|
||||||
|
if (argc != 7) {
|
||||||
|
printf("Wrong number of arguments.\n");
|
||||||
|
printf("arg1: adjacency file\n");
|
||||||
|
printf("arg2: volume file\n");
|
||||||
|
printf("arg3: output file containing particles in each zone\n");
|
||||||
|
printf("arg4: output file containing zones in each void\n");
|
||||||
|
printf("arg5: output text file\n");
|
||||||
|
printf("arg6: Density threshold (0 for no threshold)\n\n");
|
||||||
|
exit(0);
|
||||||
|
}
|
||||||
|
adjfile = argv[1];
|
||||||
|
volfile = argv[2];
|
||||||
|
zonfile = argv[3];
|
||||||
|
zonfile2 = argv[4];
|
||||||
|
txtfile = argv[5];
|
||||||
|
if (sscanf(argv[6],"%f",&voltol) == 0) {
|
||||||
|
printf("Bad density threshold.\n");
|
||||||
|
exit(0);
|
||||||
|
}
|
||||||
|
if (voltol <= 0.) {
|
||||||
|
printf("Proceeding without a density threshold.\n");
|
||||||
|
voltol = 1e30;
|
||||||
|
}
|
||||||
|
|
||||||
|
adj = fopen(adjfile, "r");
|
||||||
|
if (adj == NULL) {
|
||||||
|
printf("Unable to open %s\n",adjfile);
|
||||||
|
exit(0);
|
||||||
|
}
|
||||||
|
fread(&np,1, sizeof(int),adj);
|
||||||
|
|
||||||
|
printf("adj: %d particles\n", np);
|
||||||
|
FF;
|
||||||
|
|
||||||
|
p = (PARTICLE *)malloc(np*sizeof(PARTICLE));
|
||||||
|
/* Adjacencies*/
|
||||||
|
for (i=0;i<np;i++) {
|
||||||
|
fread(&p[i].nadj,1,sizeof(int),adj);
|
||||||
|
/* The number of adjacencies per particle */
|
||||||
|
if (p[i].nadj > 0)
|
||||||
|
p[i].adj = (int *)malloc(p[i].nadj*sizeof(int));
|
||||||
|
else p[i].adj = 0;
|
||||||
|
p[i].ncnt = 0; /* Temporarily, it's an adj counter */
|
||||||
|
}
|
||||||
|
for (i=0;i<np;i++) {
|
||||||
|
fread(&nin,1,sizeof(int),adj);
|
||||||
|
if (nin > 0)
|
||||||
|
for (k=0;k<nin;k++) {
|
||||||
|
fread(&j,1,sizeof(int),adj);
|
||||||
|
if (j < np) {
|
||||||
|
/* Set both halves of the pair */
|
||||||
|
assert(i < j);
|
||||||
|
if (p[i].ncnt == p[i].nadj)
|
||||||
|
{
|
||||||
|
p[i].adj = (int *)realloc(p[i].adj, (p[i].nadj+1)*sizeof(int));
|
||||||
|
p[i].nadj++;
|
||||||
|
}
|
||||||
|
if (p[j].ncnt == p[j].nadj)
|
||||||
|
{
|
||||||
|
p[j].adj = (int *)realloc(p[j].adj, (p[j].nadj+1)*sizeof(int));
|
||||||
|
p[j].nadj++;
|
||||||
|
}
|
||||||
|
p[i].adj[p[i].ncnt] = j;
|
||||||
|
p[j].adj[p[j].ncnt] = i;
|
||||||
|
p[i].ncnt++; p[j].ncnt++;
|
||||||
|
} else {
|
||||||
|
printf("%d: adj = %d\n",i,j);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
fclose(adj);
|
||||||
|
|
||||||
|
/* Check that we got all the pairs */
|
||||||
|
// adj = fopen(adjfile, "r");
|
||||||
|
// fread(&np,1, sizeof(int),adj);
|
||||||
|
for (i=0;i<np;i++) {
|
||||||
|
// fread(&nin,1,sizeof(int),adj); /* actually nadj */
|
||||||
|
if (p[i].ncnt != p[i].nadj) {
|
||||||
|
p[i].nadj = p[i].ncnt;
|
||||||
|
printf("We didn't get all of %d's adj's; %d != %d.\n",i,nin,p[i].nadj);
|
||||||
|
/*exit(0);*/
|
||||||
|
}
|
||||||
|
}
|
||||||
|
// fclose(adj);
|
||||||
|
|
||||||
|
/* Volumes */
|
||||||
|
vol = fopen(volfile, "r");
|
||||||
|
if (vol == NULL) {
|
||||||
|
printf("Unable to open volume file %s.\n\n",volfile);
|
||||||
|
exit(0);
|
||||||
|
}
|
||||||
|
fread(&np2,1, sizeof(int),adj);
|
||||||
|
if (np != np2) {
|
||||||
|
printf("Number of particles doesn't match! %d != %d\n",np,np2);
|
||||||
|
exit(0);
|
||||||
|
}
|
||||||
|
printf("%d particles\n", np);
|
||||||
|
FF;
|
||||||
|
for (i=0;i<np;i++) {
|
||||||
|
fread(&p[i].dens,1,sizeof(float),vol);
|
||||||
|
if ((p[i].dens < 1e-30) || (p[i].dens > 1e30)) {
|
||||||
|
printf("Whacked-out volume found, of particle %d: %f\n",i,p[i].dens);
|
||||||
|
p[i].dens = 1.;
|
||||||
|
}
|
||||||
|
p[i].dens = 1./p[i].dens; /* Get density from volume */
|
||||||
|
}
|
||||||
|
fclose(vol);
|
||||||
|
|
||||||
|
jumped = (int *)malloc(np*sizeof(int));
|
||||||
|
jumper = (int *)malloc(np*sizeof(int));
|
||||||
|
numinh = (int *)malloc(np*sizeof(int));
|
||||||
|
|
||||||
|
/* find jumper */
|
||||||
|
for (i = 0; i < np; i++) {
|
||||||
|
minvol = p[i].dens; jumper[i] = -1;
|
||||||
|
for (j=0; j< p[i].nadj; j++) {
|
||||||
|
if (p[p[i].adj[j]].dens < minvol) {
|
||||||
|
jumper[i] = p[i].adj[j];
|
||||||
|
minvol = p[jumper[i]].dens;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
numinh[i] = 0;
|
||||||
|
}
|
||||||
|
|
||||||
|
printf("About to jump ...\n"); FF;
|
||||||
|
|
||||||
|
/* Jump */
|
||||||
|
for (i = 0; i < np; i++) {
|
||||||
|
jumped[i] = i;
|
||||||
|
while (jumper[jumped[i]] > -1)
|
||||||
|
jumped[i] = jumper[jumped[i]];
|
||||||
|
numinh[jumped[i]]++;
|
||||||
|
}
|
||||||
|
printf("Post-jump ...\n"); FF;
|
||||||
|
|
||||||
|
nzones = 0;
|
||||||
|
for (i = 0; i < np; i++)
|
||||||
|
if (numinh[i] > 0) nzones++;
|
||||||
|
printf("%d initial zones found\n", nzones);
|
||||||
|
|
||||||
|
z = (ZONE *)malloc(nzones*sizeof(ZONE));
|
||||||
|
if (z == NULL) {
|
||||||
|
printf("Unable to allocate z\n");
|
||||||
|
exit(0);
|
||||||
|
}
|
||||||
|
zt = (ZONET *)malloc(nzones*sizeof(ZONET));
|
||||||
|
if (zt == NULL) {
|
||||||
|
printf("Unable to allocate zt\n");
|
||||||
|
exit(0);
|
||||||
|
}
|
||||||
|
zonenum = (int *)malloc(np*sizeof(int));
|
||||||
|
if (zonenum == NULL) {
|
||||||
|
printf("Unable to allocate zonenum\n");
|
||||||
|
exit(0);
|
||||||
|
}
|
||||||
|
|
||||||
|
h = 0;
|
||||||
|
for (i = 0; i < np; i++)
|
||||||
|
if (numinh[i] > 0) {
|
||||||
|
z[h].core = i;
|
||||||
|
zonenum[i] = h;
|
||||||
|
h++;
|
||||||
|
} else {
|
||||||
|
zonenum[i] = -1;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* Count border particles */
|
||||||
|
for (i = 0; i < np; i++)
|
||||||
|
for (j = 0; j < p[i].nadj; j++) {
|
||||||
|
testpart = p[i].adj[j];
|
||||||
|
if (jumped[i] != jumped[testpart])
|
||||||
|
zt[zonenum[jumped[i]]].nadj++;
|
||||||
|
}
|
||||||
|
|
||||||
|
for (h=0;h<nzones;h++) {
|
||||||
|
zt[h].adj = (int *)malloc(zt[h].nadj*sizeof(int));
|
||||||
|
if (zt[h].adj == NULL) {
|
||||||
|
printf("Unable to allocate %d adj's of zone %d\n",zt[h].nadj,h);
|
||||||
|
exit(0);
|
||||||
|
}
|
||||||
|
zt[h].slv = (float *)malloc(zt[h].nadj*sizeof(float));
|
||||||
|
if (zt[h].slv == NULL) {
|
||||||
|
printf("Unable to allocate %d slv's of zone %d\n",zt[h].nadj,h);
|
||||||
|
exit(0);
|
||||||
|
}
|
||||||
|
zt[h].nadj = 0;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* Find "weakest links" */
|
||||||
|
for (i = 0; i < np; i++) {
|
||||||
|
h = zonenum[jumped[i]];
|
||||||
|
for (j = 0; j < p[i].nadj; j++) {
|
||||||
|
testpart = p[i].adj[j];
|
||||||
|
if (h != zonenum[jumped[testpart]]) {
|
||||||
|
if (p[testpart].dens > p[i].dens) {
|
||||||
|
/* there could be a weakest link through testpart */
|
||||||
|
already = 0;
|
||||||
|
for (za = 0; za < zt[h].nadj; za++)
|
||||||
|
if (zt[h].adj[za] == zonenum[jumped[testpart]]) {
|
||||||
|
already = 1;
|
||||||
|
if (p[testpart].dens < zt[h].slv[za]) {
|
||||||
|
zt[h].slv[za] = p[testpart].dens;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
if (already == 0) {
|
||||||
|
zt[h].adj[zt[h].nadj] = zonenum[jumped[testpart]];
|
||||||
|
zt[h].slv[zt[h].nadj] = p[testpart].dens;
|
||||||
|
zt[h].nadj++;
|
||||||
|
}
|
||||||
|
} else { /* There could be a weakest link through i */
|
||||||
|
already = 0;
|
||||||
|
for (za = 0; za < zt[h].nadj; za++)
|
||||||
|
if (zt[h].adj[za] == zonenum[jumped[testpart]]) {
|
||||||
|
already = 1;
|
||||||
|
if (p[i].dens < zt[h].slv[za]) {
|
||||||
|
zt[h].slv[za] = p[i].dens;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
if (already == 0) {
|
||||||
|
zt[h].adj[zt[h].nadj] = zonenum[jumped[testpart]];
|
||||||
|
zt[h].slv[zt[h].nadj] = p[i].dens;
|
||||||
|
zt[h].nadj++;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
printf("Found zone adjacencies\n"); FF;
|
||||||
|
|
||||||
|
/* Free particle adjacencies */
|
||||||
|
for (i=0;i<np; i++) if (p[i].adj != 0) free(p[i].adj);
|
||||||
|
|
||||||
|
/* Use z instead of zt */
|
||||||
|
for (h=0;h<nzones;h++) {
|
||||||
|
/*printf("%d ",zt[h].nadj);*/
|
||||||
|
z[h].nadj = zt[h].nadj;
|
||||||
|
z[h].adj = (int *)malloc(zt[h].nadj*sizeof(int));
|
||||||
|
z[h].slv = (float *)malloc(zt[h].nadj*sizeof(float));
|
||||||
|
for (za = 0; za<zt[h].nadj; za++) {
|
||||||
|
z[h].adj[za] = zt[h].adj[za];
|
||||||
|
z[h].slv[za] = zt[h].slv[za];
|
||||||
|
}
|
||||||
|
free(zt[h].adj);
|
||||||
|
free(zt[h].slv);
|
||||||
|
z[h].np = numinh[z[h].core];
|
||||||
|
}
|
||||||
|
free(zt);
|
||||||
|
free(numinh);
|
||||||
|
|
||||||
|
m = (int **)malloc(nzones*sizeof(int *));
|
||||||
|
/* Not in the zone struct since it'll be freed up (contiguously, we hope)
|
||||||
|
soon */
|
||||||
|
nm = (int *)malloc(nzones*sizeof(int));
|
||||||
|
for (h=0; h<nzones; h++) {
|
||||||
|
m[h] = (int *)malloc(z[h].np*sizeof(int));
|
||||||
|
nm[h] = 0;
|
||||||
|
z[h].vol = 0.;
|
||||||
|
}
|
||||||
|
for (i=0; i<np; i++) {
|
||||||
|
h = zonenum[jumped[i]];
|
||||||
|
if (i == z[h].core) {
|
||||||
|
m[h][nm[h]] = m[h][0];
|
||||||
|
m[h][0] = i; /* Put the core particle at the top of the list */
|
||||||
|
} else {
|
||||||
|
m[h][nm[h]] = i;
|
||||||
|
}
|
||||||
|
z[h].vol += 1.0/(double)p[i].dens;
|
||||||
|
nm[h] ++;
|
||||||
|
}
|
||||||
|
free(nm);
|
||||||
|
|
||||||
|
|
||||||
|
zon = fopen(zonfile,"w");
|
||||||
|
if (zon == NULL) {
|
||||||
|
printf("Problem opening zonefile %s.\n\n",zonfile);
|
||||||
|
exit(0);
|
||||||
|
}
|
||||||
|
fwrite(&np,1,4,zon);
|
||||||
|
fwrite(&nzones,1,4,zon);
|
||||||
|
for (h=0; h<nzones; h++) {
|
||||||
|
fwrite(&(z[h].np),1,4,zon);
|
||||||
|
fwrite(m[h],z[h].np,4,zon);
|
||||||
|
free(m[h]);
|
||||||
|
}
|
||||||
|
free(m);
|
||||||
|
close(zon);
|
||||||
|
|
||||||
|
inyet = (char *)malloc(nzones*sizeof(char));
|
||||||
|
inyet2 = (char *)malloc(nzones*sizeof(char));
|
||||||
|
zonelist = (int *)malloc(nzones*sizeof(int));
|
||||||
|
zonelist2 = (int *)malloc(nzones*sizeof(int));
|
||||||
|
sorter = (double *)malloc((nzones+1)*sizeof(double));
|
||||||
|
|
||||||
|
for (h = 0; h< nzones; h++) {
|
||||||
|
inyet[h] = 0;
|
||||||
|
inyet2[h] = 0;
|
||||||
|
}
|
||||||
|
|
||||||
|
nhl = 0;
|
||||||
|
|
||||||
|
maxvol = 0.;
|
||||||
|
minvol = BIGFLT;
|
||||||
|
maxdenscontrast = 0.;
|
||||||
|
for(i=0;i<np; i++){
|
||||||
|
if (p[i].dens > maxvol) maxvol = p[i].dens;
|
||||||
|
if (p[i].dens < minvol) minvol = p[i].dens;
|
||||||
|
}
|
||||||
|
printf("Densities range from %e to %e.\n",minvol,maxvol);FF;
|
||||||
|
|
||||||
|
zon2 = fopen(zonfile2,"w");
|
||||||
|
if (zon2 == NULL) {
|
||||||
|
printf("Problem opening zonefile %s.\n\n",zonfile2);
|
||||||
|
exit(0);
|
||||||
|
}
|
||||||
|
fwrite(&nzones,1,4,zon2);
|
||||||
|
|
||||||
|
for (h = 0; h<nzones; h++) {
|
||||||
|
nhlcount = 0;
|
||||||
|
for (hl = 0; hl < nhl; hl++)
|
||||||
|
inyet[zonelist[hl]] = 0;
|
||||||
|
|
||||||
|
zonelist[0] = h;
|
||||||
|
inyet[h] = 1;
|
||||||
|
nhl = 1;
|
||||||
|
z[h].npjoin = z[h].np;
|
||||||
|
do {
|
||||||
|
/* Find the lowest-volume (highest-density) adjacency */
|
||||||
|
lowvol = BIGFLT; nl = 0; beaten = 0;
|
||||||
|
for (hl = 0; hl < nhl; hl++) {
|
||||||
|
h2 = zonelist[hl];
|
||||||
|
if (inyet[h2] == 1) { /* If it's not already identified as
|
||||||
|
an interior zone, with inyet=2 */
|
||||||
|
interior = 1;
|
||||||
|
for (za = 0; za < z[h2].nadj; za ++) {
|
||||||
|
if (inyet[z[h2].adj[za]] == 0) {
|
||||||
|
interior = 0;
|
||||||
|
if (z[h2].slv[za] == lowvol) {
|
||||||
|
link[nl] = z[h2].adj[za];
|
||||||
|
nl ++;
|
||||||
|
if (nl == NLINKS) {
|
||||||
|
printf("Too many links with the same rho_sl! Increase NLINKS from %d\n",nl);
|
||||||
|
exit(0);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
if (z[h2].slv[za] < lowvol) {
|
||||||
|
lowvol = z[h2].slv[za];
|
||||||
|
link[0] = z[h2].adj[za];
|
||||||
|
nl = 1;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
if (interior == 1) inyet[h2] = 2; /* No bordering exter. zones */
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
if (nl == 0) {
|
||||||
|
beaten = 1;
|
||||||
|
z[h].leak = maxvol;
|
||||||
|
continue;
|
||||||
|
}
|
||||||
|
|
||||||
|
if (lowvol > voltol) {
|
||||||
|
beaten = 1;
|
||||||
|
z[h].leak = lowvol;
|
||||||
|
continue;
|
||||||
|
}
|
||||||
|
|
||||||
|
for (l=0; l < nl; l++)
|
||||||
|
if (p[z[link[l]].core].dens < p[z[h].core].dens)
|
||||||
|
beaten = 1;
|
||||||
|
if (beaten == 1) {
|
||||||
|
z[h].leak = lowvol;
|
||||||
|
continue;
|
||||||
|
}
|
||||||
|
/* Add everything linked to the link(s) */
|
||||||
|
nhl2 = 0;
|
||||||
|
for (l=0; l < nl; l++) {
|
||||||
|
if (inyet2[link[l]] == 0) {
|
||||||
|
zonelist2[nhl2] = link[l];
|
||||||
|
inyet2[link[l]] = 1;
|
||||||
|
nhl2 ++;
|
||||||
|
added = 1;
|
||||||
|
while ((added == 1) && (beaten == 0)) {
|
||||||
|
added = 0;
|
||||||
|
for (hl = 0; (hl < nhl2) && (beaten == 0); hl++) {
|
||||||
|
h2 = zonelist2[hl];
|
||||||
|
if (inyet2[h2] == 1) {
|
||||||
|
interior = 1; /* Guilty until proven innocent */
|
||||||
|
for (za = 0; za < z[h2].nadj; za ++) {
|
||||||
|
link2 = z[h2].adj[za];
|
||||||
|
if ((inyet[link2]+inyet2[link2]) == 0) {
|
||||||
|
interior = 0;
|
||||||
|
if (z[h2].slv[za] <= lowvol) {
|
||||||
|
if (p[z[link2].core].dens < p[z[h].core].dens) {
|
||||||
|
beaten = 1;
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
zonelist2[nhl2] = link2;
|
||||||
|
inyet2[link2] = 1;
|
||||||
|
nhl2++;
|
||||||
|
added = 1;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
if (interior == 1) inyet2[h2] = 2;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
for (hl = 0; hl < nhl2; hl++)
|
||||||
|
inyet2[zonelist2[hl]] = 0;
|
||||||
|
|
||||||
|
/* See if there's a beater */
|
||||||
|
if (beaten == 1) {
|
||||||
|
z[h].leak = lowvol;
|
||||||
|
} else {
|
||||||
|
for (h2 = 0; h2 < nhl2; h2++) {
|
||||||
|
zonelist[nhl] = zonelist2[h2];
|
||||||
|
inyet[zonelist2[h2]] = 1;
|
||||||
|
nhl++;
|
||||||
|
z[h].npjoin += z[zonelist2[h2]].np;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
if (nhl/10000 > nhlcount) {
|
||||||
|
nhlcount = nhl/10000;
|
||||||
|
printf(" %d",nhl); FF;
|
||||||
|
}
|
||||||
|
} while((lowvol < BIGFLT) && (beaten == 0));
|
||||||
|
|
||||||
|
z[h].denscontrast = z[h].leak/p[z[h].core].dens;
|
||||||
|
if (z[h].denscontrast < 1.) z[h].denscontrast = 1.;
|
||||||
|
|
||||||
|
/* find biggest denscontrast */
|
||||||
|
if (z[h].denscontrast > maxdenscontrast) {
|
||||||
|
maxdenscontrast = (double)z[h].denscontrast;
|
||||||
|
}
|
||||||
|
|
||||||
|
/* Don't sort; want the core zone to be first */
|
||||||
|
|
||||||
|
if (nhlcount > 0) { /* Outputs the number of zones in large voids */
|
||||||
|
printf(" h%d:%d\n",h,nhl);
|
||||||
|
FF;
|
||||||
|
}
|
||||||
|
/* Calculate volume */
|
||||||
|
z[h].voljoin = 0.;
|
||||||
|
for (q = 0; q<nhl; q++) {
|
||||||
|
z[h].voljoin += z[zonelist[q]].vol;
|
||||||
|
}
|
||||||
|
|
||||||
|
z[h].nhl = nhl;
|
||||||
|
|
||||||
|
fwrite(&nhl,1,4,zon2);
|
||||||
|
fwrite(zonelist,nhl,4,zon2);
|
||||||
|
}
|
||||||
|
fclose(zon2);
|
||||||
|
|
||||||
|
printf("Maxdenscontrast = %f.\n",maxdenscontrast);
|
||||||
|
|
||||||
|
/* Assign sorter by probability (could use volume instead) */
|
||||||
|
for (h=0; h< nzones; h++)
|
||||||
|
sorter[h] = (double)z[h].denscontrast;
|
||||||
|
|
||||||
|
/* Text output file */
|
||||||
|
|
||||||
|
printf("about to sort ...\n");FF;
|
||||||
|
|
||||||
|
iord = (int *)malloc(nzones*sizeof(int));
|
||||||
|
|
||||||
|
findrtop(sorter, nzones, iord, nzones);
|
||||||
|
|
||||||
|
txt = fopen(txtfile,"w");
|
||||||
|
fprintf(txt,"%d particles, %d voids.\n", np, nzones);
|
||||||
|
fprintf(txt,"Void# FileVoid# CoreParticle CoreDens ZoneVol Zone#Part Void#Zones VoidVol Void#Part VoidDensContrast VoidProb\n");
|
||||||
|
for (h=0; h<nzones; h++) {
|
||||||
|
i = iord[h];
|
||||||
|
prob = exp(-5.12*(z[i].denscontrast-1.) - 0.8*pow(z[i].denscontrast-1.,2.8));
|
||||||
|
fprintf(txt,"%d %d %d %e %e %d %d %e %d %f %6.2e\n",
|
||||||
|
h+1, i, z[i].core, p[z[i].core].dens, z[i].vol, z[i].np, z[i].nhl, z[i].voljoin, z[i].npjoin, z[i].denscontrast, prob);
|
||||||
|
|
||||||
|
} /* h+1 to start from 1, not zero */
|
||||||
|
fclose(txt);
|
||||||
|
|
||||||
|
return(0);
|
||||||
|
}
|
||||||
47
mytools/Makefile
Normal file
47
mytools/Makefile
Normal file
|
|
@ -0,0 +1,47 @@
|
||||||
|
PROGS= showZobov
|
||||||
|
|
||||||
|
CXX=g++
|
||||||
|
CC=gcc
|
||||||
|
|
||||||
|
include config.mk
|
||||||
|
|
||||||
|
SOURCES= showZobov.cpp loadZobov.cpp zobovConf.c
|
||||||
|
|
||||||
|
LIBS= $(LDFLAGS)
|
||||||
|
|
||||||
|
all: $(PROGS)
|
||||||
|
|
||||||
|
showZobov: showZobov.o loadZobov.o zobovConf.o
|
||||||
|
|
||||||
|
depend: $(SOURCES)
|
||||||
|
@echo "[DEPENDS] $^"
|
||||||
|
@$(CC) $(CPPFLAGS) -M -MM $^ > .mydepends
|
||||||
|
|
||||||
|
distclean: clean
|
||||||
|
@rm -f .mydepends
|
||||||
|
|
||||||
|
clean:
|
||||||
|
@rm -f *.o
|
||||||
|
@rm -f $(PROGS)
|
||||||
|
|
||||||
|
.mydepends: $(SOURCES)
|
||||||
|
@touch .mydepends
|
||||||
|
@make depend
|
||||||
|
|
||||||
|
%.prog:
|
||||||
|
@echo "[L] $@"
|
||||||
|
@$(CXX) -o $@ $^ $(LIBS)
|
||||||
|
|
||||||
|
zobovConf.c zobovConf.h: showZobov.ggo Makefile
|
||||||
|
@echo "[OPT] $@"
|
||||||
|
gengetopt -i $< -f zobovConf -a zobovConf_info -F zobovConf -C
|
||||||
|
|
||||||
|
%.o: %.c
|
||||||
|
@echo "[C] $< ..."
|
||||||
|
@$(CC) -c -o $@ $< $(CFLAGS)
|
||||||
|
|
||||||
|
%.o: %.cpp
|
||||||
|
@echo "[C++] $< ..."
|
||||||
|
@$(CXX) -c -o $@ $< $(CXXFLAGS)
|
||||||
|
|
||||||
|
include .mydepends
|
||||||
6
mytools/config.mk
Normal file
6
mytools/config.mk
Normal file
|
|
@ -0,0 +1,6 @@
|
||||||
|
CC= gcc
|
||||||
|
CXX= g++
|
||||||
|
CPPFLAGS= -I$(HOME)/Science/Software/CosmoToolbox/install/include
|
||||||
|
LDFLAGS= -L$(HOME)/Science/Software/CosmoToolbox/install/lib -lCosmoTool
|
||||||
|
CXXFLAGS= $(CPPFLAGS) -ggdb -O0 -ffast-math
|
||||||
|
CFLAGS= $(CPPFLAGS) -ggdb -O0 -ffast-math
|
||||||
136
mytools/loadZobov.cpp
Normal file
136
mytools/loadZobov.cpp
Normal file
|
|
@ -0,0 +1,136 @@
|
||||||
|
#include <string>
|
||||||
|
#include <fstream>
|
||||||
|
#include <iostream>
|
||||||
|
#include <vector>
|
||||||
|
#include <cstdlib>
|
||||||
|
#include <sstream>
|
||||||
|
#include "loadZobov.hpp"
|
||||||
|
|
||||||
|
using namespace std;
|
||||||
|
|
||||||
|
bool loadZobov(const char *descName, const char *adjName, const char *volName, ZobovRep& z)
|
||||||
|
{
|
||||||
|
ifstream descFile(descName);
|
||||||
|
ifstream adjFile(adjName);
|
||||||
|
ifstream volFile(volName);
|
||||||
|
int32_t numParticles, numZones, numPinZone;
|
||||||
|
int32_t totalParticles;
|
||||||
|
int32_t numVoids;
|
||||||
|
|
||||||
|
adjFile.read((char *)&numParticles, sizeof(numParticles));
|
||||||
|
adjFile.read((char *)&numZones, sizeof(numZones));
|
||||||
|
if (!adjFile)
|
||||||
|
return false;
|
||||||
|
|
||||||
|
cout << "Number of particles = " << numParticles << endl;
|
||||||
|
cout << "Number of zones = " << numZones << endl;
|
||||||
|
|
||||||
|
totalParticles = 0;
|
||||||
|
z.allZones.resize(numZones);
|
||||||
|
for (int zone = 0; zone < numZones; zone++)
|
||||||
|
{
|
||||||
|
adjFile.read((char *)&numPinZone, sizeof(numPinZone));
|
||||||
|
if (!adjFile)
|
||||||
|
{
|
||||||
|
cout << "Problem on the zone " << zone << " / " << numZones << endl;
|
||||||
|
return false;
|
||||||
|
}
|
||||||
|
|
||||||
|
z.allZones[zone].pId.resize(numPinZone);
|
||||||
|
adjFile.read((char *)&z.allZones[zone].pId[0], sizeof(int)*numPinZone);
|
||||||
|
|
||||||
|
totalParticles += numPinZone;
|
||||||
|
}
|
||||||
|
cout << "Zoned " << totalParticles << endl;
|
||||||
|
|
||||||
|
if (totalParticles != numParticles)
|
||||||
|
{
|
||||||
|
cerr << "The numbers of particles are inconsistent ! (" << totalParticles << " vs " << numParticles << ")"<< endl;
|
||||||
|
abort();
|
||||||
|
}
|
||||||
|
|
||||||
|
volFile.read((char *)&numVoids, sizeof(numVoids));
|
||||||
|
if (!volFile)
|
||||||
|
return false;
|
||||||
|
|
||||||
|
cout << "Number of voids = " << numVoids << endl;
|
||||||
|
|
||||||
|
z.allVoids.resize(numVoids);
|
||||||
|
for (int v = 0; v < numVoids; v++)
|
||||||
|
{
|
||||||
|
int32_t numZinV;
|
||||||
|
|
||||||
|
volFile.read((char *)&numZinV, sizeof(numZinV));
|
||||||
|
if (!volFile)
|
||||||
|
return false;
|
||||||
|
|
||||||
|
z.allVoids[v].zId.resize(numZinV);
|
||||||
|
|
||||||
|
int *zId = new int[numZinV];
|
||||||
|
|
||||||
|
volFile.read((char *)zId, sizeof(int)*numZinV);
|
||||||
|
for (int k = 0; k < numZinV; k++)
|
||||||
|
z.allVoids[v].zId[k] = &z.allZones[zId[k]];
|
||||||
|
|
||||||
|
delete[] zId;
|
||||||
|
}
|
||||||
|
|
||||||
|
cout << "Loading description" << endl;
|
||||||
|
|
||||||
|
string line;
|
||||||
|
getline(descFile, line);
|
||||||
|
getline(descFile, line);
|
||||||
|
getline(descFile, line);
|
||||||
|
while (!descFile.eof())
|
||||||
|
{
|
||||||
|
istringstream lineStream(line.c_str());
|
||||||
|
int orderId, volId, coreParticle, numParticlesInZone, numZonesInVoid, numInVoid;
|
||||||
|
float coreDensity, volumeZone, volumeVoid, densityContrast;
|
||||||
|
float probability;
|
||||||
|
|
||||||
|
lineStream
|
||||||
|
>> orderId
|
||||||
|
>> volId
|
||||||
|
>> coreParticle
|
||||||
|
>> coreDensity
|
||||||
|
>> volumeZone
|
||||||
|
>> numParticlesInZone
|
||||||
|
>> numZonesInVoid
|
||||||
|
>> volumeVoid
|
||||||
|
>> numInVoid
|
||||||
|
>> densityContrast
|
||||||
|
>> probability;
|
||||||
|
if (!lineStream)
|
||||||
|
{
|
||||||
|
cerr << "Error in text stream" << endl;
|
||||||
|
abort();
|
||||||
|
}
|
||||||
|
|
||||||
|
z.allVoids[volId].proba = probability;
|
||||||
|
z.allVoids[volId].volume = volumeVoid;
|
||||||
|
z.allVoids[volId].numParticles = numInVoid;
|
||||||
|
z.allVoids[volId].coreParticle = coreParticle;
|
||||||
|
|
||||||
|
// Sanity check
|
||||||
|
int actualNumber = 0;
|
||||||
|
for (int j = 0; j < z.allVoids[volId].zId.size(); j++)
|
||||||
|
actualNumber += z.allVoids[volId].zId[j]->pId.size();
|
||||||
|
|
||||||
|
if (actualNumber != numInVoid)
|
||||||
|
{
|
||||||
|
cerr << "Sanity check failed."
|
||||||
|
<< " The number of particles in the description ("
|
||||||
|
<< numInVoid
|
||||||
|
<< ") is different from the one in the file ("
|
||||||
|
<< actualNumber << ")" << endl;
|
||||||
|
}
|
||||||
|
getline(descFile, line);
|
||||||
|
}
|
||||||
|
|
||||||
|
cout << "Done loading" << endl;
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
return true;
|
||||||
|
}
|
||||||
28
mytools/loadZobov.hpp
Normal file
28
mytools/loadZobov.hpp
Normal file
|
|
@ -0,0 +1,28 @@
|
||||||
|
#ifndef __LOAD_ZOBOV_HPP
|
||||||
|
#define __LOAD_ZOBOV_HPP
|
||||||
|
|
||||||
|
#include <vector>
|
||||||
|
|
||||||
|
struct ZobovZone
|
||||||
|
{
|
||||||
|
std::vector<int> pId;
|
||||||
|
};
|
||||||
|
|
||||||
|
struct ZobovVoid
|
||||||
|
{
|
||||||
|
std::vector<ZobovZone *> zId;
|
||||||
|
float proba;
|
||||||
|
int numParticles, coreParticle;
|
||||||
|
float volume;
|
||||||
|
};
|
||||||
|
|
||||||
|
struct ZobovRep
|
||||||
|
{
|
||||||
|
std::vector<ZobovZone> allZones;
|
||||||
|
std::vector<ZobovVoid> allVoids;
|
||||||
|
};
|
||||||
|
|
||||||
|
bool loadZobov(const char *descName,
|
||||||
|
const char *adjName, const char *volName, ZobovRep& z);
|
||||||
|
|
||||||
|
#endif
|
||||||
1068
mytools/zobovConf.c
Normal file
1068
mytools/zobovConf.c
Normal file
File diff suppressed because it is too large
Load diff
243
mytools/zobovConf.h
Normal file
243
mytools/zobovConf.h
Normal file
|
|
@ -0,0 +1,243 @@
|
||||||
|
/** @file zobovConf.h
|
||||||
|
* @brief The header file for the command line option parser
|
||||||
|
* generated by GNU Gengetopt version 2.22
|
||||||
|
* http://www.gnu.org/software/gengetopt.
|
||||||
|
* DO NOT modify this file, since it can be overwritten
|
||||||
|
* @author GNU Gengetopt by Lorenzo Bettini */
|
||||||
|
|
||||||
|
#ifndef ZOBOVCONF_H
|
||||||
|
#define ZOBOVCONF_H
|
||||||
|
|
||||||
|
/* If we use autoconf. */
|
||||||
|
#ifdef HAVE_CONFIG_H
|
||||||
|
#include "config.h"
|
||||||
|
#endif
|
||||||
|
|
||||||
|
#include <stdio.h> /* for FILE */
|
||||||
|
|
||||||
|
#ifdef __cplusplus
|
||||||
|
extern "C" {
|
||||||
|
#endif /* __cplusplus */
|
||||||
|
|
||||||
|
#ifndef ZOBOVCONF_PACKAGE
|
||||||
|
/** @brief the program name */
|
||||||
|
#define ZOBOVCONF_PACKAGE "showZobov"
|
||||||
|
#endif
|
||||||
|
|
||||||
|
#ifndef ZOBOVCONF_VERSION
|
||||||
|
/** @brief the program version */
|
||||||
|
#define ZOBOVCONF_VERSION "0"
|
||||||
|
#endif
|
||||||
|
|
||||||
|
/** @brief Where the command line options are stored */
|
||||||
|
struct zobovConf_info
|
||||||
|
{
|
||||||
|
const char *help_help; /**< @brief Print help and exit help description. */
|
||||||
|
const char *version_help; /**< @brief Print version and exit help description. */
|
||||||
|
char * desc_arg; /**< @brief The description file name for the voids. */
|
||||||
|
char * desc_orig; /**< @brief The description file name for the voids original value given at command line. */
|
||||||
|
const char *desc_help; /**< @brief The description file name for the voids help description. */
|
||||||
|
char * adj_arg; /**< @brief Adjacent file name. */
|
||||||
|
char * adj_orig; /**< @brief Adjacent file name original value given at command line. */
|
||||||
|
const char *adj_help; /**< @brief Adjacent file name help description. */
|
||||||
|
char * void_arg; /**< @brief Void/zone bind filename. */
|
||||||
|
char * void_orig; /**< @brief Void/zone bind filename original value given at command line. */
|
||||||
|
const char *void_help; /**< @brief Void/zone bind filename help description. */
|
||||||
|
int id_arg; /**< @brief Void id to extract. */
|
||||||
|
char * id_orig; /**< @brief Void id to extract original value given at command line. */
|
||||||
|
const char *id_help; /**< @brief Void id to extract help description. */
|
||||||
|
int sProba_flag; /**< @brief Sort against probability (default=off). */
|
||||||
|
const char *sProba_help; /**< @brief Sort against probability help description. */
|
||||||
|
int sVol_flag; /**< @brief Sort against volume (default=on). */
|
||||||
|
const char *sVol_help; /**< @brief Sort against volume help description. */
|
||||||
|
double mProba_arg; /**< @brief Minimal probability to accept (default='0.0'). */
|
||||||
|
char * mProba_orig; /**< @brief Minimal probability to accept original value given at command line. */
|
||||||
|
const char *mProba_help; /**< @brief Minimal probability to accept help description. */
|
||||||
|
char * ramsesDir_arg; /**< @brief Ramses base output directory. */
|
||||||
|
char * ramsesDir_orig; /**< @brief Ramses base output directory original value given at command line. */
|
||||||
|
const char *ramsesDir_help; /**< @brief Ramses base output directory help description. */
|
||||||
|
int ramsesId_arg; /**< @brief Ramses output id. */
|
||||||
|
char * ramsesId_orig; /**< @brief Ramses output id original value given at command line. */
|
||||||
|
const char *ramsesId_help; /**< @brief Ramses output id help description. */
|
||||||
|
char * configFile_arg; /**< @brief Configuration file. */
|
||||||
|
char * configFile_orig; /**< @brief Configuration file original value given at command line. */
|
||||||
|
const char *configFile_help; /**< @brief Configuration file help description. */
|
||||||
|
int quiet_flag; /**< @brief Quiet mode (default=off). */
|
||||||
|
const char *quiet_help; /**< @brief Quiet mode help description. */
|
||||||
|
char * move_arg; /**< @brief Move the center by (x,y,z) (default='(0,0,0)'). */
|
||||||
|
char * move_orig; /**< @brief Move the center by (x,y,z) original value given at command line. */
|
||||||
|
const char *move_help; /**< @brief Move the center by (x,y,z) help description. */
|
||||||
|
int galax_flag; /**< @brief Output is galax particles (default=off). */
|
||||||
|
const char *galax_help; /**< @brief Output is galax particles help description. */
|
||||||
|
char * output_arg; /**< @brief Output file (default='voidDesc'). */
|
||||||
|
char * output_orig; /**< @brief Output file original value given at command line. */
|
||||||
|
const char *output_help; /**< @brief Output file help description. */
|
||||||
|
|
||||||
|
unsigned int help_given ; /**< @brief Whether help was given. */
|
||||||
|
unsigned int version_given ; /**< @brief Whether version was given. */
|
||||||
|
unsigned int desc_given ; /**< @brief Whether desc was given. */
|
||||||
|
unsigned int adj_given ; /**< @brief Whether adj was given. */
|
||||||
|
unsigned int void_given ; /**< @brief Whether void was given. */
|
||||||
|
unsigned int id_given ; /**< @brief Whether id was given. */
|
||||||
|
unsigned int sProba_given ; /**< @brief Whether sProba was given. */
|
||||||
|
unsigned int sVol_given ; /**< @brief Whether sVol was given. */
|
||||||
|
unsigned int mProba_given ; /**< @brief Whether mProba was given. */
|
||||||
|
unsigned int ramsesDir_given ; /**< @brief Whether ramsesDir was given. */
|
||||||
|
unsigned int ramsesId_given ; /**< @brief Whether ramsesId was given. */
|
||||||
|
unsigned int configFile_given ; /**< @brief Whether configFile was given. */
|
||||||
|
unsigned int quiet_given ; /**< @brief Whether quiet was given. */
|
||||||
|
unsigned int move_given ; /**< @brief Whether move was given. */
|
||||||
|
unsigned int galax_given ; /**< @brief Whether galax was given. */
|
||||||
|
unsigned int output_given ; /**< @brief Whether output was given. */
|
||||||
|
|
||||||
|
} ;
|
||||||
|
|
||||||
|
/** @brief The additional parameters to pass to parser functions */
|
||||||
|
struct zobovConf_params
|
||||||
|
{
|
||||||
|
int override; /**< @brief whether to override possibly already present options (default 0) */
|
||||||
|
int initialize; /**< @brief whether to initialize the option structure zobovConf_info (default 1) */
|
||||||
|
int check_required; /**< @brief whether to check that all required options were provided (default 1) */
|
||||||
|
int check_ambiguity; /**< @brief whether to check for options already specified in the option structure zobovConf_info (default 0) */
|
||||||
|
int print_errors; /**< @brief whether getopt_long should print an error message for a bad option (default 1) */
|
||||||
|
} ;
|
||||||
|
|
||||||
|
/** @brief the purpose string of the program */
|
||||||
|
extern const char *zobovConf_info_purpose;
|
||||||
|
/** @brief the usage string of the program */
|
||||||
|
extern const char *zobovConf_info_usage;
|
||||||
|
/** @brief all the lines making the help output */
|
||||||
|
extern const char *zobovConf_info_help[];
|
||||||
|
|
||||||
|
/**
|
||||||
|
* The command line parser
|
||||||
|
* @param argc the number of command line options
|
||||||
|
* @param argv the command line options
|
||||||
|
* @param args_info the structure where option information will be stored
|
||||||
|
* @return 0 if everything went fine, NON 0 if an error took place
|
||||||
|
*/
|
||||||
|
int zobovConf (int argc, char * const *argv,
|
||||||
|
struct zobovConf_info *args_info);
|
||||||
|
|
||||||
|
/**
|
||||||
|
* The command line parser (version with additional parameters - deprecated)
|
||||||
|
* @param argc the number of command line options
|
||||||
|
* @param argv the command line options
|
||||||
|
* @param args_info the structure where option information will be stored
|
||||||
|
* @param override whether to override possibly already present options
|
||||||
|
* @param initialize whether to initialize the option structure my_args_info
|
||||||
|
* @param check_required whether to check that all required options were provided
|
||||||
|
* @return 0 if everything went fine, NON 0 if an error took place
|
||||||
|
* @deprecated use zobovConf_ext() instead
|
||||||
|
*/
|
||||||
|
int zobovConf2 (int argc, char * const *argv,
|
||||||
|
struct zobovConf_info *args_info,
|
||||||
|
int override, int initialize, int check_required);
|
||||||
|
|
||||||
|
/**
|
||||||
|
* The command line parser (version with additional parameters)
|
||||||
|
* @param argc the number of command line options
|
||||||
|
* @param argv the command line options
|
||||||
|
* @param args_info the structure where option information will be stored
|
||||||
|
* @param params additional parameters for the parser
|
||||||
|
* @return 0 if everything went fine, NON 0 if an error took place
|
||||||
|
*/
|
||||||
|
int zobovConf_ext (int argc, char * const *argv,
|
||||||
|
struct zobovConf_info *args_info,
|
||||||
|
struct zobovConf_params *params);
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Save the contents of the option struct into an already open FILE stream.
|
||||||
|
* @param outfile the stream where to dump options
|
||||||
|
* @param args_info the option struct to dump
|
||||||
|
* @return 0 if everything went fine, NON 0 if an error took place
|
||||||
|
*/
|
||||||
|
int zobovConf_dump(FILE *outfile,
|
||||||
|
struct zobovConf_info *args_info);
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Save the contents of the option struct into a (text) file.
|
||||||
|
* This file can be read by the config file parser (if generated by gengetopt)
|
||||||
|
* @param filename the file where to save
|
||||||
|
* @param args_info the option struct to save
|
||||||
|
* @return 0 if everything went fine, NON 0 if an error took place
|
||||||
|
*/
|
||||||
|
int zobovConf_file_save(const char *filename,
|
||||||
|
struct zobovConf_info *args_info);
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Print the help
|
||||||
|
*/
|
||||||
|
void zobovConf_print_help(void);
|
||||||
|
/**
|
||||||
|
* Print the version
|
||||||
|
*/
|
||||||
|
void zobovConf_print_version(void);
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Initializes all the fields a zobovConf_params structure
|
||||||
|
* to their default values
|
||||||
|
* @param params the structure to initialize
|
||||||
|
*/
|
||||||
|
void zobovConf_params_init(struct zobovConf_params *params);
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Allocates dynamically a zobovConf_params structure and initializes
|
||||||
|
* all its fields to their default values
|
||||||
|
* @return the created and initialized zobovConf_params structure
|
||||||
|
*/
|
||||||
|
struct zobovConf_params *zobovConf_params_create(void);
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Initializes the passed zobovConf_info structure's fields
|
||||||
|
* (also set default values for options that have a default)
|
||||||
|
* @param args_info the structure to initialize
|
||||||
|
*/
|
||||||
|
void zobovConf_init (struct zobovConf_info *args_info);
|
||||||
|
/**
|
||||||
|
* Deallocates the string fields of the zobovConf_info structure
|
||||||
|
* (but does not deallocate the structure itself)
|
||||||
|
* @param args_info the structure to deallocate
|
||||||
|
*/
|
||||||
|
void zobovConf_free (struct zobovConf_info *args_info);
|
||||||
|
|
||||||
|
/**
|
||||||
|
* The config file parser (deprecated version)
|
||||||
|
* @param filename the name of the config file
|
||||||
|
* @param args_info the structure where option information will be stored
|
||||||
|
* @param override whether to override possibly already present options
|
||||||
|
* @param initialize whether to initialize the option structure my_args_info
|
||||||
|
* @param check_required whether to check that all required options were provided
|
||||||
|
* @return 0 if everything went fine, NON 0 if an error took place
|
||||||
|
* @deprecated use zobovConf_config_file() instead
|
||||||
|
*/
|
||||||
|
int zobovConf_configfile (char * const filename,
|
||||||
|
struct zobovConf_info *args_info,
|
||||||
|
int override, int initialize, int check_required);
|
||||||
|
|
||||||
|
/**
|
||||||
|
* The config file parser
|
||||||
|
* @param filename the name of the config file
|
||||||
|
* @param args_info the structure where option information will be stored
|
||||||
|
* @param params additional parameters for the parser
|
||||||
|
* @return 0 if everything went fine, NON 0 if an error took place
|
||||||
|
*/
|
||||||
|
int zobovConf_config_file (char * const filename,
|
||||||
|
struct zobovConf_info *args_info,
|
||||||
|
struct zobovConf_params *params);
|
||||||
|
|
||||||
|
/**
|
||||||
|
* Checks that all the required options were specified
|
||||||
|
* @param args_info the structure to check
|
||||||
|
* @param prog_name the name of the program that will be used to print
|
||||||
|
* possible errors
|
||||||
|
* @return
|
||||||
|
*/
|
||||||
|
int zobovConf_required (struct zobovConf_info *args_info,
|
||||||
|
const char *prog_name);
|
||||||
|
|
||||||
|
|
||||||
|
#ifdef __cplusplus
|
||||||
|
}
|
||||||
|
#endif /* __cplusplus */
|
||||||
|
#endif /* ZOBOVCONF_H */
|
||||||
121
readfiles.c
Normal file
121
readfiles.c
Normal file
|
|
@ -0,0 +1,121 @@
|
||||||
|
#include <stdio.h>
|
||||||
|
#include <stdlib.h>
|
||||||
|
|
||||||
|
#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; i<np; i++) {
|
||||||
|
(*p)[i] = (float *)malloc(3*sizeof(float));
|
||||||
|
if ((*p)[i] == NULL) {
|
||||||
|
printf("Unable to allocate particle array in readfiles!\n");
|
||||||
|
fflush(stdout);
|
||||||
|
exit(0);
|
||||||
|
}
|
||||||
|
(*p)[i][0] = ptemp[i];
|
||||||
|
}
|
||||||
|
fread(&dum,1,4,pos);
|
||||||
|
fread(&dum,1,4,pos);
|
||||||
|
fread(ptemp,np,4,pos);
|
||||||
|
for (i=0; i<np; i++) (*p)[i][1] = ptemp[i];
|
||||||
|
fread(&dum,1,4,pos);
|
||||||
|
fread(&dum,1,4,pos);
|
||||||
|
fread(ptemp,np,4,pos);
|
||||||
|
for (i=0; i<np; i++) (*p)[i][2] = ptemp[i];
|
||||||
|
fread(&dum,1,4,pos);
|
||||||
|
|
||||||
|
fclose(pos);
|
||||||
|
free(ptemp);
|
||||||
|
|
||||||
|
/* Get into physical units (Mpc/h) */
|
||||||
|
|
||||||
|
|
||||||
|
printf("%f\n",fact);fflush(stdout);
|
||||||
|
for (i=0; i<np; i++) DL (*p)[i][d] *= fact;
|
||||||
|
|
||||||
|
/* Test range -- can comment out */
|
||||||
|
xmin = BF; xmax = -BF; ymin = BF; ymax = -BF; zmin = BF; zmax = -BF;
|
||||||
|
for (i=0; i<np;i++) {
|
||||||
|
if ((*p)[i][0]<xmin) xmin = (*p)[i][0]; if ((*p)[i][0]>xmax) xmax = (*p)[i][0];
|
||||||
|
if ((*p)[i][1]<ymin) ymin = (*p)[i][1]; if ((*p)[i][1]>ymax) ymax = (*p)[i][1];
|
||||||
|
if ((*p)[i][2]<zmin) zmin = (*p)[i][2]; 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<np; i++) DL (*v)[d][i] *= fact;
|
||||||
|
|
||||||
|
/* Test range -- can comment out */
|
||||||
|
xmin = BF; xmax = -BF; ymin = BF; ymax = -BF; zmin = BF; zmax = -BF;
|
||||||
|
for (i=0; i<np;i++) {
|
||||||
|
if ((*v)[0][i] < xmin) xmin = (*v)[0][i]; if ((*v)[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);
|
||||||
|
}
|
||||||
25
readme.txt
Normal file
25
readme.txt
Normal file
|
|
@ -0,0 +1,25 @@
|
||||||
|
This is the readme file for ZOBOV (ZOnes Bordering On Voidness),
|
||||||
|
version 1.0 (Jan 5, 2008) a void-finding algorithm for sets of
|
||||||
|
particles (e.g. a cosmological N-body simulation) by Mark Neyrinck.
|
||||||
|
It is an inversion of the halo-finding algorithm VOBOZ (VOronoi BOund
|
||||||
|
Zones), which was also written by Mark Neyrinck, with the guidance of
|
||||||
|
Nick Gnedin and Andrew Hamilton.
|
||||||
|
|
||||||
|
The gz file which once contained this file should include the
|
||||||
|
following files:
|
||||||
|
|
||||||
|
Makefile jozov.c readme.txt voz1b1.c voztie.c
|
||||||
|
findrtop.c readfiles.c voz.h vozinit.c vozutil.c
|
||||||
|
|
||||||
|
For instructions on installing and running ZOBOV, please see
|
||||||
|
http://www.ifa.hawaii.edu/~neyrinck/zobov/. To use ZOBOV, you will
|
||||||
|
need to download the Qhull package, at http://www.qhull.org. For help
|
||||||
|
in understanding ZOBOV, please see the paper describing it
|
||||||
|
(arXiv:0712.3049), If these resources are insufficient, feel free to
|
||||||
|
email Mark.Neyrinck@colorado.edu with questions.
|
||||||
|
|
||||||
|
This is free software. It may be freely copied, modified, and
|
||||||
|
redistributed, as long as ZOBOV and its author are acknowledged.
|
||||||
|
There is no warranty or other guarantee of fitness for ZOBOV; it is
|
||||||
|
provided "as is." I welcome bug reports (but do not guarantee their
|
||||||
|
fixing), which should be sent to Mark.Neyrinck@colorado.edu.
|
||||||
8
voz.h
Normal file
8
voz.h
Normal file
|
|
@ -0,0 +1,8 @@
|
||||||
|
#define MAXVERVER 2000
|
||||||
|
#define NGUARD 42 /*Actually, the number of SPACES between guard points
|
||||||
|
in each dim */
|
||||||
|
|
||||||
|
typedef struct Partadj {
|
||||||
|
int nadj;
|
||||||
|
int *adj;
|
||||||
|
} PARTADJ;
|
||||||
334
voz1b1.c
Normal file
334
voz1b1.c
Normal file
|
|
@ -0,0 +1,334 @@
|
||||||
|
#include "qhull_a.h"
|
||||||
|
#include "voz.h"
|
||||||
|
|
||||||
|
#define DL for (d=0;d<3;d++)
|
||||||
|
#define BF 1e30
|
||||||
|
|
||||||
|
int delaunadj (coordT *points, int nvp, int nvpbuf, int nvpall, PARTADJ **adjs);
|
||||||
|
int vorvol (coordT *deladjs, coordT *points, pointT *intpoints, int numpoints, float *vol);
|
||||||
|
int posread(char *posfile, float ***p, float fact);
|
||||||
|
|
||||||
|
int main(int argc, char *argv[]) {
|
||||||
|
int exitcode;
|
||||||
|
int i, j, np;
|
||||||
|
float **r;
|
||||||
|
coordT rtemp[3], *parts;
|
||||||
|
coordT deladjs[3*MAXVERVER], points[3*MAXVERVER];
|
||||||
|
pointT intpoints[3*MAXVERVER];
|
||||||
|
FILE *pos, *out;
|
||||||
|
char *posfile, outfile[80], *suffix;
|
||||||
|
PARTADJ *adjs;
|
||||||
|
float *vols;
|
||||||
|
float predict, xmin,xmax,ymin,ymax,zmin,zmax;
|
||||||
|
int *orig;
|
||||||
|
|
||||||
|
int isitinbuf;
|
||||||
|
char isitinmain, d;
|
||||||
|
int numdiv, nvp, nvpall, nvpbuf;
|
||||||
|
float width, width2, totwidth, totwidth2, bf, s, g;
|
||||||
|
float border, boxsize;
|
||||||
|
float c[3];
|
||||||
|
int b[3];
|
||||||
|
double totalvol;
|
||||||
|
|
||||||
|
if (argc != 9) {
|
||||||
|
printf("Wrong number of arguments.\n");
|
||||||
|
printf("arg1: position file\n");
|
||||||
|
printf("arg2: border size\n");
|
||||||
|
printf("arg3: box size\n");
|
||||||
|
printf("arg4: suffix\n");
|
||||||
|
printf("arg5: number of divisions\n");
|
||||||
|
printf("arg6-8: b[0-2]\n\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);
|
||||||
|
}
|
||||||
|
|
||||||
|
/* Boxsize should be the range in r, yielding a range 0-1 */
|
||||||
|
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; i<np;i++) {
|
||||||
|
if (r[i][0]<xmin) xmin = r[i][0]; if (r[i][0]>xmax) xmax = r[i][0];
|
||||||
|
if (r[i][1]<ymin) ymin = r[i][1]; if (r[i][1]>ymax) ymax = r[i][1];
|
||||||
|
if (r[i][2]<zmin) zmin = r[i][2]; if (r[i][2]>zmax) zmax = r[i][2];
|
||||||
|
}
|
||||||
|
printf("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("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(0);
|
||||||
|
}
|
||||||
|
|
||||||
|
DL c[d] = ((float)b[d]+0.5)*width;
|
||||||
|
printf("c: %f,%f,%f\n",c[0],c[1],c[2]);
|
||||||
|
/* Assign temporary array*/
|
||||||
|
nvpbuf = 0; /* Number of particles to tesselate, including
|
||||||
|
buffer */
|
||||||
|
nvp = 0; /* Without the buffer */
|
||||||
|
for (i=0; i<np; i++) {
|
||||||
|
isitinbuf = 1;
|
||||||
|
isitinmain = 1;
|
||||||
|
DL {
|
||||||
|
rtemp[d] = (double)r[i][d] - (double)c[d];
|
||||||
|
if (rtemp[d] > 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 = (coordT *)malloc(3*nvpbuf*sizeof(coordT));
|
||||||
|
orig = (int *)malloc(nvpbuf*sizeof(int));
|
||||||
|
|
||||||
|
if (parts == NULL) {
|
||||||
|
printf("Unable to allocate parts\n");
|
||||||
|
fflush(stdout);
|
||||||
|
}
|
||||||
|
if (orig == NULL) {
|
||||||
|
printf("Unable to allocate orig\n");
|
||||||
|
fflush(stdout);
|
||||||
|
}
|
||||||
|
|
||||||
|
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<np; i++) {
|
||||||
|
isitinmain = 1;
|
||||||
|
DL {
|
||||||
|
rtemp[d] = r[i][d] - c[d];
|
||||||
|
if (rtemp[d] > 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("nvp = %d\n",nvp);
|
||||||
|
printf("x: %f,%f; y: %f,%f; z:%f,%f\n",xmin,xmax,ymin,ymax,zmin,zmax);
|
||||||
|
nvpbuf = nvp;
|
||||||
|
for (i=0; i<np; i++) {
|
||||||
|
isitinbuf = 1;
|
||||||
|
DL {
|
||||||
|
rtemp[d] = r[i][d] - c[d];
|
||||||
|
if (rtemp[d] > 0.5) rtemp[d] --;
|
||||||
|
if (rtemp[d] < -0.5) rtemp[d] ++;
|
||||||
|
isitinbuf = isitinbuf && (fabs(rtemp[d])<totwidth2);
|
||||||
|
}
|
||||||
|
if ((isitinbuf > 0) &&
|
||||||
|
((fabs(rtemp[0])>width2)||(fabs(rtemp[1])>width2)||(fabs(rtemp[2])>width2))) {
|
||||||
|
|
||||||
|
/*printf("%3.3f ",sqrt(rtemp[0]*rtemp[0] + rtemp[1]*rtemp[1] +
|
||||||
|
rtemp[2]*rtemp[2]));
|
||||||
|
printf("|%2.2f,%2.2f,%2.2f,%f,%f",r[i][0],r[i][1],r[i][2],width2,totwidth2);*/
|
||||||
|
parts[3*nvpbuf] = rtemp[0];
|
||||||
|
parts[3*nvpbuf+1] = rtemp[1];
|
||||||
|
parts[3*nvpbuf+2] = rtemp[2];
|
||||||
|
orig[nvpbuf] = i;
|
||||||
|
|
||||||
|
nvpbuf++;
|
||||||
|
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("nvpbuf = %d\n",nvpbuf);
|
||||||
|
printf("x: %f,%f; y: %f,%f; z:%f,%f\n",xmin,xmax,ymin,ymax,zmin,zmax);
|
||||||
|
nvpall = nvpbuf;
|
||||||
|
predict = pow(width+2.*bf,3)*(float)np;
|
||||||
|
printf("There should be ~ %f points; there are %d\n",predict,nvpbuf);
|
||||||
|
|
||||||
|
for (i=0;i<np;i++) free(r[i]);
|
||||||
|
free(r);
|
||||||
|
|
||||||
|
/* Add guard points */
|
||||||
|
for (i=0; i<NGUARD+1; i++) {
|
||||||
|
for (j=0; j<NGUARD+1; j++) {
|
||||||
|
/* Bottom */
|
||||||
|
parts[3*nvpall] = -width2 + (float)i * s;
|
||||||
|
parts[3*nvpall+1] = -width2 + (float)j * s;
|
||||||
|
parts[3*nvpall+2] = -width2 - g;
|
||||||
|
nvpall++;
|
||||||
|
/* Top */
|
||||||
|
parts[3*nvpall] = -width2 + (float)i * s;
|
||||||
|
parts[3*nvpall+1] = -width2 + (float)j * s;
|
||||||
|
parts[3*nvpall+2] = width2 + g;
|
||||||
|
nvpall++;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
for (i=0; i<NGUARD+1; i++) { /* Don't want to overdo the corners*/
|
||||||
|
for (j=0; j<NGUARD+1; j++) {
|
||||||
|
parts[3*nvpall] = -width2 + (float)i * s;
|
||||||
|
parts[3*nvpall+1] = -width2 - g;
|
||||||
|
parts[3*nvpall+2] = -width2 + (float)j * s;
|
||||||
|
nvpall++;
|
||||||
|
|
||||||
|
parts[3*nvpall] = -width2 + (float)i * s;
|
||||||
|
parts[3*nvpall+1] = width2 + g;
|
||||||
|
parts[3*nvpall+2] = -width2 + (float)j * s;
|
||||||
|
nvpall++;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
for (i=0; i<NGUARD+1; i++) {
|
||||||
|
for (j=0; j<NGUARD+1; j++) {
|
||||||
|
parts[3*nvpall] = -width2 - g;
|
||||||
|
parts[3*nvpall+1] = -width2 + (float)i * s;
|
||||||
|
parts[3*nvpall+2] = -width2 + (float)j * s;
|
||||||
|
nvpall++;
|
||||||
|
|
||||||
|
parts[3*nvpall] = width2 + g;
|
||||||
|
parts[3*nvpall+1] = -width2 + (float)i * s;
|
||||||
|
parts[3*nvpall+2] = -width2 + (float)j * s;
|
||||||
|
nvpall++;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
xmin = BF; xmax = -BF; ymin = BF; ymax = -BF; zmin = BF; zmax = -BF;
|
||||||
|
for (i=nvpbuf;i<nvpall;i++) {
|
||||||
|
if (parts[3*i] < xmin) xmin = parts[3*i];
|
||||||
|
if (parts[3*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("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);
|
||||||
|
exitcode = delaunadj(parts, nvp, nvpbuf, nvpall, &adjs);
|
||||||
|
|
||||||
|
/* Now calculate volumes*/
|
||||||
|
printf("Now finding volumes ...\n"); fflush(stdout);
|
||||||
|
vols = (float *)malloc(nvp*sizeof(float));
|
||||||
|
|
||||||
|
for (i=0; i<nvp; i++) { /* Just the original particles
|
||||||
|
Assign adjacency coordinate array*/
|
||||||
|
/* Volumes */
|
||||||
|
for (j = 0; j < adjs[i].nadj; j++)
|
||||||
|
DL {
|
||||||
|
deladjs[3*j + d] = parts[3*adjs[i].adj[j]+d] - parts[3*i+d];
|
||||||
|
if (deladjs[3*j+d] < -0.5) deladjs[3*j+d]++;
|
||||||
|
if (deladjs[3*j+d] > 0.5) deladjs[3*j+d]--;
|
||||||
|
}
|
||||||
|
|
||||||
|
exitcode = vorvol(deladjs, points, intpoints, adjs[i].nadj, &(vols[i]));
|
||||||
|
vols[i] *= (float)np;
|
||||||
|
if (i % 1000 == 0)
|
||||||
|
printf("%d: %d, %f\n",i,adjs[i].nadj,vols[i]);
|
||||||
|
}
|
||||||
|
|
||||||
|
/* Get the adjacencies back to their original values */
|
||||||
|
|
||||||
|
for (i=0; i<nvp; i++)
|
||||||
|
for (j = 0; j < adjs[i].nadj; j++)
|
||||||
|
adjs[i].adj[j] = orig[adjs[i].adj[j]];
|
||||||
|
|
||||||
|
totalvol = 0.;
|
||||||
|
for (i=0;i<nvp; i++) {
|
||||||
|
totalvol += (double)vols[i];
|
||||||
|
}
|
||||||
|
printf("Average volume = %g\n",totalvol/(float)nvp);
|
||||||
|
|
||||||
|
/* Now the output!
|
||||||
|
First number of particles */
|
||||||
|
sprintf(outfile,"part.%s.%02d.%02d.%02d",suffix,b[0],b[1],b[2]);
|
||||||
|
|
||||||
|
printf("Output to %s\n\n",outfile);
|
||||||
|
out = fopen(outfile,"w");
|
||||||
|
if (out == NULL) {
|
||||||
|
printf("Unable to open %s\n",outfile);
|
||||||
|
exit(0);
|
||||||
|
}
|
||||||
|
fwrite(&np,1, sizeof(int),out);
|
||||||
|
fwrite(&nvp,1, sizeof(int),out);
|
||||||
|
printf("nvp = %d\n",nvp);
|
||||||
|
|
||||||
|
/* Tell us where the original particles were */
|
||||||
|
fwrite(orig,sizeof(int),nvp,out);
|
||||||
|
/* Volumes*/
|
||||||
|
fwrite(vols,sizeof(float),nvp,out);
|
||||||
|
/* Adjacencies */
|
||||||
|
for (i=0;i<nvp;i++) {
|
||||||
|
fwrite(&(adjs[i].nadj),1,sizeof(int),out);
|
||||||
|
if (adjs[i].nadj > 0)
|
||||||
|
fwrite(adjs[i].adj,adjs[i].nadj,sizeof(int),out);
|
||||||
|
else printf("0");
|
||||||
|
}
|
||||||
|
fclose(out);
|
||||||
|
|
||||||
|
return(0);
|
||||||
|
}
|
||||||
153
vozinit.c
Normal file
153
vozinit.c
Normal file
|
|
@ -0,0 +1,153 @@
|
||||||
|
#include "qhull_a.h"
|
||||||
|
#include "voz.h"
|
||||||
|
|
||||||
|
#define DL for (d=0;d<3;d++)
|
||||||
|
#define BF 1e30
|
||||||
|
|
||||||
|
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[80], systemstr[90], *suffix;
|
||||||
|
float xmin,xmax,ymin,ymax,zmin,zmax;
|
||||||
|
|
||||||
|
int isitinbuf;
|
||||||
|
char isitinmain, d;
|
||||||
|
int numdiv;
|
||||||
|
int nvp, nvpall, nvpbuf, nvpmin, nvpmax, nvpbufmin, nvpbufmax; /* yes, the insurance */
|
||||||
|
float width, width2, totwidth, totwidth2, bf, s, g;
|
||||||
|
float border, boxsize;
|
||||||
|
float c[3];
|
||||||
|
int b[3];
|
||||||
|
|
||||||
|
if (argc != 6) {
|
||||||
|
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\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 < 2) {
|
||||||
|
printf("Cannot have a number of divisions less than 2. Resetting to 2:\n");
|
||||||
|
numdiv = 2;
|
||||||
|
}
|
||||||
|
|
||||||
|
suffix = argv[5];
|
||||||
|
|
||||||
|
/* Read the position file */
|
||||||
|
np = posread(posfile,&rfloat,1./boxsize);
|
||||||
|
/* Boxsize should be the range in r, yielding a range 0-1 */
|
||||||
|
|
||||||
|
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;
|
||||||
|
|
||||||
|
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<np; i++) {
|
||||||
|
isitinbuf = 1; isitinmain = 1;
|
||||||
|
for (d=0; d<3; d++) {
|
||||||
|
rtemp[d] = rfloat[i][d] - c[d];
|
||||||
|
if (rtemp[d] > 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("b=(%d,%d,%d), c=(%f,%f,%f), nvp=%d, nvpbuf=%d\n",
|
||||||
|
b[0],b[1],b[2],c[0],c[1],c[2],nvp,nvpbuf);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
printf("Nvp range: %d,%d\n",nvpmin,nvpmax);
|
||||||
|
printf("Nvpbuf range: %d,%d\n",nvpbufmin,nvpbufmax);
|
||||||
|
|
||||||
|
/* 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(0);
|
||||||
|
}
|
||||||
|
fprintf(scr,"#!/bin/bash -f\n");
|
||||||
|
for (b[0]=0;b[0]<numdiv; b[0]++) {
|
||||||
|
for (b[1] = 0; b[1] < numdiv; b[1]++) {
|
||||||
|
for (b[2] = 0; b[2] < numdiv; b[2]++) {
|
||||||
|
fprintf(scr,"./voz1b1 %s %f %f %s %d %d %d %d\n",
|
||||||
|
posfile,border,boxsize,suffix,numdiv,b[0],b[1],b[2]);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
fprintf(scr,"voztie %d %s\n",numdiv,suffix);
|
||||||
|
fclose(scr);
|
||||||
|
|
||||||
|
sprintf(systemstr,"chmod u+x %s",scrfile);
|
||||||
|
system(systemstr);
|
||||||
|
|
||||||
|
return(0);
|
||||||
|
}
|
||||||
182
voztie.c
Normal file
182
voztie.c
Normal file
|
|
@ -0,0 +1,182 @@
|
||||||
|
#include <math.h>
|
||||||
|
#include <stdio.h>
|
||||||
|
#include <stdlib.h>
|
||||||
|
#include "voz.h"
|
||||||
|
|
||||||
|
int main(int argc, char *argv[]) {
|
||||||
|
|
||||||
|
FILE *part, *adj, *vol;
|
||||||
|
char partfile[80], *suffix, adjfile[80], volfile[80];
|
||||||
|
float *vols, volstemp;
|
||||||
|
|
||||||
|
PARTADJ *adjs;
|
||||||
|
|
||||||
|
int numdiv,np,np2,na;
|
||||||
|
|
||||||
|
int i,j,k,p,nout;
|
||||||
|
int nvp,npnotdone,nvpmax, nvpsum, *orig;
|
||||||
|
double avgnadj, avgvol;
|
||||||
|
|
||||||
|
if (argc != 3) {
|
||||||
|
printf("Wrong number of arguments.\n");
|
||||||
|
printf("arg1: number of divisions (default 2)\n");
|
||||||
|
printf("arg2: suffix describing this run\n\n");
|
||||||
|
exit(0);
|
||||||
|
}
|
||||||
|
if (sscanf(argv[1],"%d",&numdiv) != 1) {
|
||||||
|
printf("That's no number of divisions; try again.\n");
|
||||||
|
exit(0);
|
||||||
|
}
|
||||||
|
if (numdiv < 2) {
|
||||||
|
printf("Cannot have a number of divisions less than 2. Resetting to 2:\n");
|
||||||
|
numdiv = 2;
|
||||||
|
}
|
||||||
|
suffix = argv[2];
|
||||||
|
|
||||||
|
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,"part.%s.%02d.%02d.%02d",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);
|
||||||
|
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);
|
||||||
|
|
||||||
|
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 = (int *)malloc(nvpmax*sizeof(int));
|
||||||
|
if (orig == NULL) printf("Couldn't allocate orig.\n");
|
||||||
|
if ((vols == NULL) || (orig == NULL) || (adjs == NULL)) {
|
||||||
|
printf("Not enough memory to allocate. Exiting.\n");
|
||||||
|
exit(0);
|
||||||
|
}
|
||||||
|
for (p=0;p<np;p++)
|
||||||
|
vols[p] = -1.;
|
||||||
|
|
||||||
|
for (i = 0; i < numdiv; i++) {
|
||||||
|
for (j = 0; j < numdiv; j++) {
|
||||||
|
for (k = 0; k < numdiv; k++) {
|
||||||
|
sprintf(partfile,"part.%s.%02d.%02d.%02d",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);
|
||||||
|
/*printf("nvp = %d\n",nvp);fflush(stdout);*/
|
||||||
|
|
||||||
|
nvpsum += nvp;
|
||||||
|
|
||||||
|
fread(orig,nvp,sizeof(int),part);
|
||||||
|
for (p=0;p<nvp;p++) {
|
||||||
|
fread(&volstemp,1,sizeof(float),part);
|
||||||
|
if (vols[orig[p]] > -1.)
|
||||||
|
if (fabs(vols[orig[p]]-volstemp)/volstemp > 1e-3) {
|
||||||
|
printf("Inconsistent volumes for p. %d: (%g,%g)!\n",
|
||||||
|
orig[p],vols[orig[p]],volstemp);
|
||||||
|
exit(0);
|
||||||
|
}
|
||||||
|
vols[orig[p]] = volstemp;
|
||||||
|
}
|
||||||
|
|
||||||
|
for (p=0;p<nvp;p++) {
|
||||||
|
fread(&na,1,sizeof(int),part);
|
||||||
|
if (na > 0) {
|
||||||
|
adjs[orig[p]].nadj = na;
|
||||||
|
adjs[orig[p]].adj = (int *)malloc(na*sizeof(int));
|
||||||
|
if (adjs[orig[p]].adj == NULL) {
|
||||||
|
printf("Couldn't allocate adjs[orig[%d]].adj.\n",p);
|
||||||
|
exit(0);
|
||||||
|
}
|
||||||
|
fread(adjs[orig[p]].adj,na,sizeof(int),part);
|
||||||
|
} else {
|
||||||
|
printf("0"); fflush(stdout);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
fclose(part);
|
||||||
|
printf("%d ",k);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
printf("\n");
|
||||||
|
npnotdone = 0; avgnadj = 0.; avgvol = 0.;
|
||||||
|
for (p=0;p<np;p++) {
|
||||||
|
if (vols[p] == -1.) npnotdone++;
|
||||||
|
avgnadj += (double)(adjs[p].nadj);
|
||||||
|
avgvol += (double)(vols[p]);
|
||||||
|
}
|
||||||
|
if (npnotdone > 0)
|
||||||
|
printf("%d particles not done!\n");
|
||||||
|
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,"adj%s.dat",suffix);
|
||||||
|
sprintf(volfile,"vol%s.dat",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 */
|
||||||
|
for (i=0;i<np;i++)
|
||||||
|
fwrite(&adjs[i].nadj,1,sizeof(int),adj);
|
||||||
|
|
||||||
|
/* Now the lists of adjacencies (without double counting) */
|
||||||
|
for (i=0;i<np;i++)
|
||||||
|
if (adjs[i].nadj > 0) {
|
||||||
|
nout = 0;
|
||||||
|
for (j=0;j<adjs[i].nadj; j++) if (adjs[i].adj[j] > i) nout++;
|
||||||
|
fwrite(&nout,1,sizeof(int),adj);
|
||||||
|
for (j=0;j<adjs[i].nadj; j++)
|
||||||
|
if (adjs[i].adj[j] > i)
|
||||||
|
fwrite(&(adjs[i].adj[j]),1,sizeof(int),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);
|
||||||
|
|
||||||
|
return(0);
|
||||||
|
}
|
||||||
241
vozutil.c
Normal file
241
vozutil.c
Normal file
|
|
@ -0,0 +1,241 @@
|
||||||
|
#include "qhull_a.h"
|
||||||
|
#include "voz.h"
|
||||||
|
|
||||||
|
#define FOREACHvertex2_(vertices) FOREACHsetelement_(vertexT, vertices2,vertex2)
|
||||||
|
|
||||||
|
/*char qh_version[] = "user_eg 3.1 2001/10/04"; */
|
||||||
|
|
||||||
|
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 (coordT *points, int nvp, int nvpbuf, int nvpall, PARTADJ **adjs) {
|
||||||
|
|
||||||
|
int dim= 3; /* dimension of points */
|
||||||
|
boolT 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;
|
||||||
|
|
||||||
|
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");
|
||||||
|
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, NULL, &numfacets, &numsimplicial,
|
||||||
|
&totneighbors, &numridges, &numcoplanars, &numtricoplanars);
|
||||||
|
qh_vertexneighbors();
|
||||||
|
vertices= qh_facetvertices (qh facet_list, NULL, NULL);
|
||||||
|
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<adjst.nadj; i++)
|
||||||
|
if (adjst.adj[i] != adjst.adj[i-1]) {
|
||||||
|
if (adjst.adj[i] >= nvpbuf) {
|
||||||
|
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<adjst.nadj; i++)
|
||||||
|
if (adjst.adj[i] != adjst.adj[i-1]) {
|
||||||
|
(*adjs)[ver].adj[count] = adjst.adj[i];
|
||||||
|
count++;
|
||||||
|
}
|
||||||
|
(*adjs)[ver].nadj = count;
|
||||||
|
} else {
|
||||||
|
printf("Number of adjacencies %d < 4, particle %d -> %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 (coordT *deladjs, coordT *points, pointT *intpoints, int numpoints, float *vol) {
|
||||||
|
int dim= 3; /* dimension of points */
|
||||||
|
boolT 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 */
|
||||||
|
|
||||||
|
coordT *point, *normp, *coordp, *feasiblep, *deladj;
|
||||||
|
int i, j, k;
|
||||||
|
boolT zerodiv;
|
||||||
|
float runsum;
|
||||||
|
char region;
|
||||||
|
/*coordT *points;
|
||||||
|
pointT *intpoints;*/
|
||||||
|
|
||||||
|
/* make point array from adjacency coordinates (add offset)*/
|
||||||
|
/*points = (coordT *)malloc(4*numpoints*sizeof(coordT));
|
||||||
|
if (points == NULL) {
|
||||||
|
printf("Unable to allocate points\n");
|
||||||
|
exit(0);
|
||||||
|
}*/
|
||||||
|
for (i=0; i<numpoints; i++) {
|
||||||
|
runsum = 0.;
|
||||||
|
deladj = deladjs + 3*i;
|
||||||
|
point = points + 4*i;
|
||||||
|
for (j=0;j<3;j++) {
|
||||||
|
runsum += deladj[j]*deladj[j];
|
||||||
|
point[j] = deladj[j];
|
||||||
|
}
|
||||||
|
point[3] = -0.5*runsum;
|
||||||
|
}
|
||||||
|
sprintf (flags, "qhull H0");
|
||||||
|
|
||||||
|
exitcode= qh_new_qhull (4, numpoints, points, ismalloc,
|
||||||
|
flags, outfile, errfile);
|
||||||
|
|
||||||
|
numpoints = 0;
|
||||||
|
if (!exitcode) { /* if no error */
|
||||||
|
FORALLfacets {
|
||||||
|
numpoints++;
|
||||||
|
}
|
||||||
|
/* Now we know how many points */
|
||||||
|
/*intpoints = (pointT *)malloc(dim*numpoints*sizeof(pointT));
|
||||||
|
if (intpoints == NULL) {
|
||||||
|
printf("Unable to allocate intpoints\n");
|
||||||
|
exit(0);
|
||||||
|
}*/
|
||||||
|
|
||||||
|
j = 0;
|
||||||
|
FORALLfacets {
|
||||||
|
if (!qh feasible_point) {
|
||||||
|
fprintf (stdout, "qhull input error (qh_printafacet): option 'Fp' needs qh feasible_point\n");
|
||||||
|
qh_errexit( qh_ERRinput, NULL, NULL);
|
||||||
|
}
|
||||||
|
point= coordp= intpoints + j*3;
|
||||||
|
j++;
|
||||||
|
normp= facet->normal;
|
||||||
|
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;
|
||||||
|
}
|
||||||
Loading…
Add table
Add a link
Reference in a new issue