From 5ba4109f0021e7b2c9c66821461742a339e07355 Mon Sep 17 00:00:00 2001 From: pants Date: Wed, 7 Dec 2016 13:29:51 -0500 Subject: added support for hyperuniform lattices using existing hard-sphere jamming routine --- .gitignore | 2 + makefile | 32 +- makefile_ept | 31 -- makefile_hal | 31 -- src/fracture.c | 6 +- src/fracture.h | 5 +- src/graph_create.c | 2 + src/graph_genfunc.c | 53 +- src/spheres_poly/box.cpp | 1157 +++++++++++++++++++++++++++++++++++++++++ src/spheres_poly/box.h | 169 ++++++ src/spheres_poly/event.cpp | 65 +++ src/spheres_poly/event.h | 58 +++ src/spheres_poly/grid_field.h | 148 ++++++ src/spheres_poly/heap.cpp | 186 +++++++ src/spheres_poly/heap.h | 42 ++ src/spheres_poly/neighbor.cpp | 40 ++ src/spheres_poly/sphere.cpp | 66 +++ src/spheres_poly/sphere.h | 44 ++ src/spheres_poly/spheres.cpp | 64 +++ src/spheres_poly/spheres.h | 4 + src/spheres_poly/vector.h | 471 +++++++++++++++++ 21 files changed, 2551 insertions(+), 125 deletions(-) delete mode 100644 makefile_ept delete mode 100644 makefile_hal create mode 100644 src/spheres_poly/box.cpp create mode 100644 src/spheres_poly/box.h create mode 100644 src/spheres_poly/event.cpp create mode 100644 src/spheres_poly/event.h create mode 100644 src/spheres_poly/grid_field.h create mode 100644 src/spheres_poly/heap.cpp create mode 100644 src/spheres_poly/heap.h create mode 100644 src/spheres_poly/neighbor.cpp create mode 100644 src/spheres_poly/sphere.cpp create mode 100644 src/spheres_poly/sphere.h create mode 100644 src/spheres_poly/spheres.cpp create mode 100644 src/spheres_poly/spheres.h create mode 100644 src/spheres_poly/vector.h diff --git a/.gitignore b/.gitignore index 2eae92a..bebedfe 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,5 @@ bin/* obj/* data/raw/* +data/raw +makefile_* diff --git a/makefile b/makefile index d788743..ed1142d 100644 --- a/makefile +++ b/makefile @@ -1,12 +1,18 @@ CC = clang +CCP = clang++ CFLAGS = -g -Os -O3 -Wall -fno-strict-aliasing -Wstrict-overflow -Wno-missing-field-initializers -fPIC -flto -fopenmp -march=native +CPFLAGS = -g -Os -O3 -Wall -fno-strict-aliasing -Wstrict-overflow -Wno-missing-field-initializers -fPIC -flto -fopenmp -march=native +CCPFLAGS = -g -Os -O3 -Wall -fno-strict-aliasing -Wstrict-overflow -Wno-missing-field-initializers -fPIC -flto -shared -fopenmp -march=native LDFLAGS = -lc -lcblas -llapack -ldl -lpthread -lcholmod -lamd -lcolamd -lsuitesparseconfig -lcamd -lccolamd -lm -lrt -lmetis -lgsl -lprofiler -ltcmalloc OBJ = data bound_set correlations factor_update graph_genfunc net net_voltages net_currents net_conductivity net_fracture get_dual_clusters break_edge graph_components gen_laplacian geometry gen_voltcurmat graph_create graph_free fortune/edgelist fortune/geometry fortune/heap fortune/main fortune/output fortune/voronoi fortune/memory rand -BIN = corr_test fracture anal_process big_anal_process cen_anal_process +OBJP = spheres_poly/box spheres_poly/event spheres_poly/heap spheres_poly/neighbor spheres_poly/sphere +OBJCCP = spheres_poly/spheres +BIN = corr_test fracture anal_process big_anal_process cen_anal_process long_anal_process -all: opt ${OBJ:%=obj/%.o} ${BIN:%=obj/%.o} ${BIN:%=bin/%} + +all: opt ${OBJ:%=obj/%.o} ${OBJP:%=obj/%.o} ${BIN:%=obj/%.o} ${BIN:%=bin/%} opt: @echo build options: @@ -14,18 +20,26 @@ opt: @echo "CFLAGS = ${CFLAGS}" @echo "LDFLAGS = ${LDFLAGS}" -obj/%.o: src/%.c +${OBJ:%=obj/%.o} ${BIN:%=obj/%.o}: obj/%.o: src/%.c @echo CC -c -o $@ - @${CC} -c -o $@ $< ${CFLAGS} + @${CC} -c ${CFLAGS} $< -o $@ + +${OBJP:%=obj/%.o}: obj/%.o: src/%.cpp + @echo CCP -c -o $@ + @${CCP} -c ${CPFLAGS} $< -o $@ + +${OBJCCP:%=obj/%.o}: obj/%.o: src/%.cpp ${OBJP:%=obj/%.o} + @echo ${CCP} ${CCPFLAGS} ${OBJP:%=obj/%.o} $< -o $@ + @${CCP} ${CCPFLAGS} ${OBJP:%=obj/%.o} -fuse-ld=gold $< -o $@ -bin/%: obj/%.o ${OBJ:%=obj/%.o} - @echo CC -o $@ - @${CC} -o $@ $< ${OBJ:%=obj/%.o} -fuse-ld=gold ${CFLAGS} ${LDFLAGS} +bin/%: obj/%.o ${OBJ:%=obj/%.o} ${OBJCCP:%=obj/%.o} + @echo ${CC} ${OBJ:%=obj/%.o} ${OBJCCP:%=obj/%.o} -fuse-ld=gold ${CFLAGS} ${LDFLAGS} $< -o $@ + @${CC} ${OBJ:%=obj/%.o} ${OBJCCP:%=obj/%.o} -fuse-ld=gold ${CFLAGS} ${LDFLAGS} $< -o $@ clean: @echo cleaning: - @echo rm -f ${OBJ:%=obj/%} ${BIN:%=obj/%.o} ${BIN:%=bin/%} - @rm -f ${OBJ:%=obj/%.o} ${BIN:%=obj/%.o} ${BIN:%=bin/%} + @echo rm -f ${OBJ:%=obj/%.o} ${OBJP:%=obj/%.o} ${OBJCCP:%=obj/%.o} ${BIN:%=obj/%.o} ${BIN:%=bin/%} + @rm -f ${OBJ:%=obj/%.o} ${OBJP:%=obj/%.o} ${OBJCCP:%=obj/%.o} ${BIN:%=obj/%.o} ${BIN:%=bin/%} .PHONY: all clean diff --git a/makefile_ept b/makefile_ept deleted file mode 100644 index c38af25..0000000 --- a/makefile_ept +++ /dev/null @@ -1,31 +0,0 @@ - -CC = clang -CFLAGS = -g -Os -O3 -Wall -fno-strict-aliasing -Wstrict-overflow -Wno-missing-field-initializers -fPIC -flto -fopenmp=libiomp5 -march=native -I/usr/include/suitesparse -I/usr/lib/gcc/x86_64-linux-gnu/5/include -LDFLAGS = -lc -lcblas -llapack -ldl -lpthread -lcholmod -lamd -lcolamd -lsuitesparseconfig -lcamd -lccolamd -lm -lrt -lmetis -lgsl -lprofiler -ltcmalloc - -OBJ = data bound_set correlations factor_update graph_genfunc net net_voltages net_currents net_conductivity net_fracture get_dual_clusters break_edge graph_components gen_laplacian geometry gen_voltcurmat graph_create graph_free fortune/edgelist fortune/geometry fortune/heap fortune/main fortune/output fortune/voronoi fortune/memory rand -BIN = corr_test fracture anal_process big_anal_process cen_anal_process long_anal_process - -all: opt ${OBJ:%=obj/%.o} ${BIN:%=obj/%.o} ${BIN:%=bin/%} - -opt: - @echo build options: - @echo "CC = ${CC}" - @echo "CFLAGS = ${CFLAGS}" - @echo "LDFLAGS = ${LDFLAGS}" - -obj/%.o: src/%.c - @echo CC -c -o $@ - @${CC} -c -o $@ $< ${CFLAGS} - -bin/%: obj/%.o ${OBJ:%=obj/%.o} - @echo CC -o $@ - @${CC} -o $@ $< ${OBJ:%=obj/%.o} -fuse-ld=gold ${CFLAGS} ${LDFLAGS} - -clean: - @echo cleaning: - @echo rm -f ${OBJ:%=obj/%} ${BIN:%=obj/%.o} ${BIN:%=bin/%} - @rm -f ${OBJ:%=obj/%.o} ${BIN:%=obj/%.o} ${BIN:%=bin/%} - -.PHONY: all clean - diff --git a/makefile_hal b/makefile_hal deleted file mode 100644 index e12d1d7..0000000 --- a/makefile_hal +++ /dev/null @@ -1,31 +0,0 @@ - -CC = clang -CFLAGS = -g -Os -O3 -Wall -fno-strict-aliasing -Wstrict-overflow -Wno-missing-field-initializers -fPIC -flto -fopenmp -LDFLAGS = -lc -lcblas -llapack -ldl -lpthread -lcholmod -lamd -lcolamd -lsuitesparseconfig -lcamd -lccolamd -lm -lrt -lmetis -lgsl -lprofiler #-ltcmalloc - -OBJ = data bound_set correlations factor_update graph_genfunc net net_voltages net_currents net_conductivity net_fracture get_dual_clusters break_edge graph_components gen_laplacian geometry gen_voltcurmat graph_create graph_free fortune/edgelist fortune/geometry fortune/heap fortune/main fortune/output fortune/voronoi fortune/memory rand -BIN = corr_test fracture anal_process big_anal_process cen_anal_process long_anal_process - -all: opt ${OBJ:%=obj/%.o} ${BIN:%=obj/%.o} ${BIN:%=bin/%} - -opt: - @echo build options: - @echo "CC = ${CC}" - @echo "CFLAGS = ${CFLAGS}" - @echo "LDFLAGS = ${LDFLAGS}" - -obj/%.o: src/%.c - @echo CC -c -o $@ - @${CC} -c -o $@ $< ${CFLAGS} - -bin/%: obj/%.o ${OBJ:%=obj/%.o} - @echo CC -o $@ - @${CC} -o $@ $< ${OBJ:%=obj/%.o} -fuse-ld=gold ${CFLAGS} ${LDFLAGS} - -clean: - @echo cleaning: - @echo rm -f ${OBJ:%=obj/%} ${BIN:%=obj/%.o} ${BIN:%=bin/%} - @rm -f ${OBJ:%=obj/%.o} ${BIN:%=obj/%.o} ${BIN:%=bin/%} - -.PHONY: all clean - diff --git a/src/fracture.c b/src/fracture.c index d155b02..bb7701b 100644 --- a/src/fracture.c +++ b/src/fracture.c @@ -105,8 +105,12 @@ int main(int argc, char *argv[]) { lattice = SQUARE_LATTICE; lattice_c = 's'; break; + case 2: + lattice = VORONOI_HYPERUNIFORM_LATTICE; + lattice_c = 'h'; + break; default: - printf("lattice specifier must be 0 (VORONOI_LATTICE) or 1 (SQUARE_LATTICE).\n"); + printf("lattice specifier must be 0 (VORONOI_LATTICE), 1 (SQUARE_LATTICE), or 2 (VORONOI_HYPERUNIFORM_LATTICE).\n"); exit(EXIT_FAILURE); } break; diff --git a/src/fracture.h b/src/fracture.h index f4f7334..8f883bc 100644 --- a/src/fracture.h +++ b/src/fracture.h @@ -33,7 +33,8 @@ typedef enum lattice_t { VORONOI_LATTICE, - SQUARE_LATTICE + SQUARE_LATTICE, + VORONOI_HYPERUNIFORM_LATTICE } lattice_t; typedef enum bound_t { @@ -185,3 +186,5 @@ bool set_connected(const cholmod_sparse *laplacian, uint_t *marks, int vertex, i unsigned long int rand_seed(); long double rand_dist_pow(const gsl_rng *r, double beta); + +double *spheres(int N, double bidispersityratio, double bidispersityfraction, double maxpf, double maxpressure, double maxcollisionrate); diff --git a/src/graph_create.c b/src/graph_create.c index 635b12b..b071c49 100644 --- a/src/graph_create.c +++ b/src/graph_create.c @@ -672,6 +672,8 @@ graph_t *graph_create(lattice_t lattice, bound_t bound, uint_t L, bool dual, cho if (bound == EMBEDDED_BOUND) side_bounds = true; else side_bounds = false; return ini_square_network(L, bound, side_bounds, c); + case (VORONOI_HYPERUNIFORM_LATTICE): + return ini_voro_graph(L, bound, dual, genfunc_hyperuniform, c); default: printf("lattice type unsupported\n"); exit(EXIT_FAILURE); diff --git a/src/graph_genfunc.c b/src/graph_genfunc.c index 54d9daa..8dcd0df 100644 --- a/src/graph_genfunc.c +++ b/src/graph_genfunc.c @@ -20,58 +20,7 @@ double g(double rho, double dist) { double *genfunc_hyperuniform(unsigned int L, bound_t boundary, gsl_rng *r, unsigned int *num) { *num = 2 * pow(L / 2, 2); - // necessary to prevent crashing when underflow occurs - gsl_set_error_handler_off(); - - double *lattice = (double *)malloc(2 * (*num) * sizeof(double)); - double rho = *num; - unsigned int to_gen = *num; - unsigned int start = 0; - - if (boundary == EMBEDDED_BOUND) { - for (unsigned int i = 0; i < L / 2; i++) { - lattice[2 * i + 1] = 0; - lattice[2 * i] = (2. * i + 1.) / L; - - lattice[L + 2 * i + 1] = 1; - lattice[L + 2 * i] = (2. * i + 1.) / L; - - lattice[2 * L + 2 * i + 1] = (2. * i + 1.) / L; - lattice[2 * L + 2 * i] = 0; - - lattice[3 * L + 2 * i + 1] = (2. * i + 1.) / L; - lattice[3 * L + 2 * i] = 1; - } - - to_gen -= 2 * L; - start = 2 * L; - } - - for (unsigned int i = 0; i < to_gen; i++) { - bool reject = true; - double x, y; - while(reject) { - x = gsl_ran_flat(r, 0, 1); - y = gsl_ran_flat(r, 0, 1); - double pp = 1; - for (unsigned int j = 0; j < i; j++) { - double ds0, ds1, ds2, ds3, ds4, ds5, ds6, ds7, ds8; - ds0 = pow(x-lattice[2*j],2)+pow(y-lattice[2*j+1],2); - ds1 = pow(x-lattice[2*j] + 1,2)+pow(y-lattice[2*j+1],2); - ds2 = pow(x-lattice[2*j] - 1,2)+pow(y-lattice[2*j+1],2); - ds3 = pow(x-lattice[2*j],2)+pow(y-lattice[2*j+1] + 1,2); - ds4 = pow(x-lattice[2*j],2)+pow(y-lattice[2*j+1] - 1,2); - ds5 = pow(x-lattice[2*j] + 1,2)+pow(y-lattice[2*j+1] + 1,2); - ds6 = pow(x-lattice[2*j] + 1,2)+pow(y-lattice[2*j+1] - 1,2); - ds7 = pow(x-lattice[2*j] - 1,2)+pow(y-lattice[2*j+1] + 1,2); - ds8 = pow(x-lattice[2*j] - 1,2)+pow(y-lattice[2*j+1] - 1,2); - pp *= g(rho, ds0) * g(rho, ds1) * g(rho, ds2) * g(rho, ds3) * g(rho, ds4) * g(rho, ds5) * g(rho, ds6) * g(rho, ds7) * g(rho, ds8); - } - if (pp > gsl_ran_flat(r, 0, 1)) reject = false; - } - lattice[2*start + 2 * i] = x; - lattice[2*start + 2 * i + 1] = y; - } + double *lattice = spheres(*num, 0.8, 0.5, 0.9, 100, 100000); return lattice; } diff --git a/src/spheres_poly/box.cpp b/src/spheres_poly/box.cpp new file mode 100644 index 0000000..716090d --- /dev/null +++ b/src/spheres_poly/box.cpp @@ -0,0 +1,1157 @@ +/* + Updated July 24, 2009 to include hardwall boundary condition option, as + well as polydisperse spheres. + + Packing of hard spheres via molecular dynamics + Developed by Monica Skoge, 2006, Princeton University + Contact: Monica Skoge (mskoge.princeton.edu) with questions + This code may be used, modified and distributed freely. + Please cite: + + "Packing Hyperspheres in High-Dimensional Euclidean Spaces" + M. Skoge, A. Donev, F. H. Stillinger and S. Torquato, 2006 + + if you use these codes. +*/ + +#include "box.h" +#include +#include +#include +#include +#include + +//============================================================== +//============================================================== +// Class Box: Fills box with hardspheres to given packing fraction +// and evolves spheres using molecular dynamics! +//============================================================== +//============================================================== + + +//============================================================== +// Constructor +//============================================================== +box::box(int N_i, double r_i, double growthrate_i, double maxpf_i, + double bidispersityratio_i, double bidispersityfraction_i, + double massratio_i, int hardwallBC_i): + r(r_i), + N(N_i), + growthrate(growthrate_i), + h(N_i+1), + maxpf(maxpf_i), + bidispersityratio(bidispersityratio_i), + bidispersityfraction(bidispersityfraction_i), + massratio(massratio_i), + hardwallBC(hardwallBC_i) +{ + ngrids = Optimalngrids2(r); + cells.set_size(ngrids); + + s = new sphere[N]; + binlist = new int[N]; + x = new vector[N]; + h.s = s; + + gtime = 0.; + rtime = 0.; + ncollisions = 0; + ntransfers = 0; + nchecks = 0; + neventstot = 0; + ncycles = 0; + xmomentum = 0.; + pressure = 0.; + collisionrate = 0.; + + cells.set_size(ngrids); + cells.initialize(-1); // initialize cells to -1 + srandom(::time(0)); // initialize the random number generator + for (int i=0; i> s[i].r; // read in radius + infile >> s[i].gr; // read in growth rate + infile >> s[i].m; // read in mass + for (int k=0; k> s[i].x[k]; // read in position + } + infile.close(); +} + + +//============================================================== +// Recreates all N spheres at random positions +//============================================================== +void box::RecreateSpheres(const char* filename, double temp) +{ + ReadPositions(filename); // reads in positions of spheres + VelocityGiver(temp); // gives spheres initial velocities + AssignCells(); // assigns spheres to cells + SetInitialEvents(); +} + + +//============================================================== +// Creates all N spheres at random positions +//============================================================== +void box::CreateSpheres(double temp) +{ + int Ncurrent = 0; + for(int i=0; i xrand; // random new position vector + double d = 0.; + double radius; + double growth_rate; + double mass; + int species; + + if (Ncurrent < bidispersityfraction*N) + { + radius = r; + growth_rate = growthrate; + mass = 1.; + species = 1; + } + else + { + radius = r*bidispersityratio; + growth_rate = growthrate*bidispersityratio; + mass = massratio; + species = 2; + } + + while (counter<1000) + { + keeper = 1; + + for(int k=0; k SIZE*SIZE/4.) + { + if (xrand[k] > SIZE/2.) // look at right image + d += (xrand[k] - (s[i].x[k]+SIZE))* + (xrand[k] - (s[i].x[k]+SIZE)); + else // look at left image + d += (xrand[k] - (s[i].x[k]-SIZE))* + (xrand[k] - (s[i].x[k]-SIZE)); + } + else + d += (xrand[k] - s[i].x[k])*(xrand[k] - s[i].x[k]); + } + } + if (d <= (radius + s[i].r)*(radius + s[i].r)) // overlapping! + { + keeper = 0; + counter++; + break; + } + } + + if (hardwallBC) + { + for (int k=0; k= 1000) + { + std::cout << "counter >= 1000" << std::endl; + exit(-1); + } + } + + // now convert xrand into index vector for cells + vector cell; + cell = vector::integer(xrand*((double)(ngrids))/SIZE); + + s[Ncurrent] = sphere(Ncurrent, xrand, cell, gtime, radius, growth_rate, + mass, species); + + //first check to see if entry at cell + if (cells.get(cell) == -1) //if yes, add Ncurrent to cells gridfield + cells.get(cell) = Ncurrent; + + else // if no, add i to right place in binlist + { + int iterater = cells.get(cell); // now iterate through to end and add Ncurrent + int pointer = iterater; + while (iterater != -1) + { + pointer = iterater; + iterater = binlist[iterater]; + } + binlist[pointer] = Ncurrent; + } + +} + + +//============================================================== +// Assign cells to spheres read in from existing configuration +//============================================================== +void box::AssignCells() +{ + for (int i=0; i cell; + cell = vector::integer(s[i].x*((double)(ngrids))/SIZE); + s[i].cell = cell; + + //first check to see if entry at cell + if (cells.get(cell) == -1) //if yes, add Ncurrent to cells gridfield + cells.get(cell) = i; + + else // if no, add i to right place in binlist + { + int iterater = cells.get(cell); // now iterate through to end and add Ncurrent + int pointer = iterater; + while (iterater != -1) + { + pointer = iterater; + iterater = binlist[iterater]; + } + binlist[pointer] = i; + } + } +} + + +//============================================================== +// Velocity Giver, assigns initial velocities from Max/Boltz dist. +//============================================================== +void box::VelocityGiver(double T) +{ + for (int i=0; i1)) + { + event o = event(outgrowtime,i,INF-1); + return o; + } + + if ((c.time < t.time)&&(c.j == INF)) // next event is check at DBL infinity + return c; + else if (c.time < t.time) // next event is collision! + { + CollisionChecker(c); + return c; + } + else // next event is transfer! + return t; +} + + +//============================================================== +// Checks events of predicted collision partner to keep collisions +// symmetric +//============================================================== +void box::CollisionChecker(event c) +{ + int i = c.i; + int j = c.j; + event cj(c.time,j,i,c.v*(-1)); + + // j should have NO event before collision with i! + if (!(c.time < s[j].nextevent.time)) + std::cout << i << " " << j << " error collchecker, s[j].nextevent.time= " << s[j].nextevent.time << " " << s[j].nextevent.j << ", c.time= " << c.time << std::endl; + + int k = s[j].nextevent.j; + if ((k < N) && (k!=i)) // j's next event was collision so give k a check + s[k].nextevent.j = INF; + + // give collision cj to j + s[j].nextevent = cj; + h.upheap(h.index[j]); +} + + +//============================================================== +// Find next transfer for sphere i +//============================================================== +event box::FindNextTransfer(int i) +{ + double ttime = dblINF; + int wallindex = INF; // -(k+1) left wall, (k+1) right wall + + vector xi = s[i].x + s[i].v*(gtime - s[i].lutime); + vector vi = s[i].v; + + for (int k=0; k0) // will hit right wall, need to check if last wall + { + if ((hardwallBC)&&(s[i].cell[k] == ngrids - 1)) + newtime = ((double)(s[i].cell[k]+1)*SIZE/((double)(ngrids)) + - (xi[k]+s[i].r+s[i].gr*gtime))/(vi[k]+s[i].gr); + else + newtime = ((double)(s[i].cell[k]+1)*SIZE/((double)(ngrids)) + - xi[k])/(vi[k]); + + if (newtime vl, vector vr, + neighbor& operation) +{ + vector cell = s[i].cell; + + // now iterate through nearest neighbors + vector offset; // nonnegative neighbor offset + vector pboffset; // nearest image offset + + vector grid; + + int ii ; + + grid=vl; + while(1) + { + //if (vr[0] > 1) + //std::cout << grid << "..." << cell+grid << "\n"; + for(int k=0; k=ngrids) // out of bounds to right + pboffset[k] = 1; + else + pboffset[k] = 0; + } + int j = cells.get((cell+offset)%ngrids); + while(j!=-1) + { + operation.Operation(j,pboffset); + j = binlist[j]; + } + + // A. Donev: + // This code makes this loop dimension-independent + // It is basically a flattened-out loop nest of depth DIM + for(ii=0;ii=DIM) break; + } +} + + +//============================================================== +// PredictCollision +//============================================================== +void box::PredictCollision(int i, int j, vector pboffset, + double& ctime, int& cpartner, + vector& cpartnerpboffset) +{ + double ctimej; + + if (i!=j) + { + ctimej = CalculateCollision(i,j,pboffset.Double())+gtime; + + if (ctimej < gtime) + std::cout << "error in find collision ctimej < 0" << std::endl; + + if ((ctimej < ctime)&&(ctimej < s[j].nextevent.time)) + { + ctime = ctimej; + cpartner = j; + cpartnerpboffset = pboffset; + } + } +} + + +//============================================================== +// Find next collision +//============================================================== +event box::FindNextCollision(int i) +{ + collision cc(i, this); + + vector vl, vr; + + for (int k=0; k pboffset) +{ + if ((hardwallBC)&&(pboffset.norm_squared() > 1E-12)) + { + return dblINF; + } + else + { + // calculate updated position and velocity of i and j + vector xi = s[i].x + s[i].v*(gtime - s[i].lutime); + vector vi = s[i].v; + vector xj = s[j].x + pboffset*SIZE + s[j].v*(gtime - s[j].lutime); + vector vj = s[j].v; + + double r_now_i = s[i].r + gtime*s[i].gr; + double r_now_j = s[j].r + gtime*s[j].gr; + + double A,B,C; + A = vector::norm_squared(vi - vj) - (s[i].gr+s[j].gr)*(s[i].gr+s[j].gr); + B = vector::dot(xi - xj, vi - vj) - (r_now_i+r_now_j)*(s[i].gr+s[j].gr); + C = vector::norm_squared(xi - xj) - (r_now_i+r_now_j)*(r_now_i+r_now_j); + + if (C < -1E-12*(r_now_i+r_now_j)) + { + std::cout << "error, " << i << " and " << j << + " are overlapping at time "<< gtime << std::cout; + std::cout << "A, B, C = " << A << " " << " " << B << + " " << " " << C << std::endl; + if (CheckSphereDiameters()>0) + std::cout << "a sphere has grown greater than unit cell" << + std::endl; + else + std::cout << "unknown error" << std::cout; + exit(-1); + } + + return QuadraticFormula(A, B, C); + } +} + + +//============================================================== +// Quadratic Formula ax^2 + bx + c = 0 +//============================================================== + double box::QuadraticFormula(double a, double b, double c) +{ + double x = dblINF; + double xpos; + double xneg; + double det = b*b - a*c; + + if (c <= 0.) + { + if(b < 0.) // spheres already overlapping and approaching + { + //std::cout << "spheres overlapping and approaching" << std::endl; + //std::cout << "# events= " << neventstot << std::endl; + x = 0.; + } + } + else if (det > -10.*DBL_EPSILON) + { + if (det < 0.) // determinant can be very small for double roots + det = 0.; + if (b < 0.) + x = c/(-b + sqrt(det)); + else if ((a < 0.)&&(b > 0.)) + x = -(b + sqrt(det))/a; + else + x = dblINF; + } + return x; +} + + +//============================================================== +// Returns first event +//============================================================== +void box::ProcessEvent() +{ + neventstot++; + // Extract first event from heap + int i = h.extractmax(); + event e = s[i].nextevent; // current event + event f; // replacement event + + if ((e.j>=0)&&(e.j vipar; // parallel comp. vi + vector vjpar; // parallel comp. vj + vector viperp; // perpendicular comp. vi + vector vjperp; // perpendicular comp. vj + + // make unit vector out of displacement vector + vector dhat; + dhat = s[i].x - s[j].x - v.Double()*SIZE; // using image of j!! + double dhatmagnitude = sqrt(dhat.norm_squared()); + dhat /= dhatmagnitude; + + vipar = dhat*vector::dot(s[i].v, dhat); + vjpar = dhat*vector::dot(s[j].v, dhat); + viperp = s[i].v - vipar; + vjperp = s[j].v - vjpar; + + s[i].v = vjpar + dhat*(s[i].gr+s[j].gr)*2 + viperp; + s[j].v = vipar - dhat*(s[i].gr+s[j].gr)*2 + vjperp; + + // momentum exchange + double xvelocity; // exchanged velocity + xvelocity = vector::dot(s[i].v - s[j].v, dhat) - (s[i].gr+s[j].gr); + xmomentum += xvelocity*dhatmagnitude*s[i].m*s[j].m*2/(s[i].m+s[j].m); +} + + +//============================================================== +// Transfer, takes care of boundary events too +//============================================================= +void box::Transfer(event e) +{ + gtime = e.time; + int i = e.i; + int j = e.j; + int k=0; // dimension perpendicular to wall it crosses + + // update position and lutime (velocity doesn't change) + s[i].x += s[i].v*(gtime-s[i].lutime); + s[i].lutime = gtime; + + vector celli; // new cell for i + celli = s[i].cell; // this is not redundant + + // update cell + if (j>N+DIM+1) // right wall + { + k = j-N-DIM-2; + celli[k] = s[i].cell[k] + 1; + + if (hardwallBC) + { + // if in right-most cell, reflect back + if (s[i].cell[k] == ngrids - 1) + s[i].v[k] = -s[i].v[k] - 4.*s[i].gr; + else + { + if (ngrids>1) + UpdateCell(i, celli); + } + } + else + { + // if in right-most cell, translate x and cell + if (s[i].cell[k] == ngrids - 1) + { + s[i].x[k] -= SIZE; + celli[k] -= ngrids; + } + if (ngrids>1) + UpdateCell(i, celli); + } + } + else if (j1) + UpdateCell(i, celli); + } + } + else + { + // if in left-most cell, translate x and cell + if (s[i].cell[k] == 0) + { + s[i].x[k] += SIZE; + celli[k] += ngrids; + } + if (ngrids>1) + UpdateCell(i, celli); + } + } + else + std::cout << "error in Transfer" << std::endl; + +} + + +//============================================================== +// Updates cell of a sphere to time +//============================================================= +void box::UpdateCell(int i, vector& celli) +{ + if (celli == s[i].cell) + std::cout << "error in update cell..shouldn't be the same" << std::endl; + + // delete i from cell array at cell + + if (cells.get(s[i].cell) == i) + { + if (binlist[i] == -1) + cells.get(s[i].cell) = -1; + else + { + cells.get(s[i].cell) = binlist[i]; + binlist[i] = -1; + } + } + + else if (cells.get(s[i].cell) == -1) + { + std::cout << "error " << i << " not in claimed cell UpdateCell" << std::endl; + OutputCells(); + } + + else // if no, find i in binlist + { + int iterater = cells.get(s[i].cell); + int pointer = iterater; + while ((iterater != i)&&(iterater != -1)) + { + pointer = iterater; + iterater = binlist[iterater]; + } + if (iterater == -1) // got to end of list without finding i + { + std::cout << "problem " << i << " wasn't in claimed, cell iterater = -1" << std::endl; + OutputCells(); + } + else // we found i! + { + binlist[pointer] = binlist[i]; + binlist[i] = -1; + } + } + + // now add i to cell array at celli + s[i].cell = celli; + + //first check to see if entry at celli + if (cells.get(celli) == -1) //if yes, add i to cells gridfield + cells.get(celli) = i; + else // if no, add i to right place in binlist + { + int iterater = cells.get(celli); // now iterate through to end and add i + int pointer = iterater; + while (iterater != -1) // find the end of the list + { + pointer = iterater; + iterater = binlist[iterater]; + } + binlist[pointer] = i; + binlist[i] = -1; // redundant + } +} + + +//============================================================== +// Output event heap...purely used for debugging +//============================================================== +void box::OutputEvents() +{ + h.print(); +} + + +//============================================================== +// Output positions of spheres and their cells...purely used for debugging +//============================================================== +void box::OutputCells() +{ + for (int i=0; i 1./ngrids){ + offender = i; + break; + } + } + return offender; +} + + +//============================================================== +// Change ngrids +//============================================================== +void box::ChangeNgrids(int newngrids) +{ + cells.set_size(newngrids); + cells.initialize(-1); // initialize cells to -1 + for (int i=0; i +#include + +#include "vector.h" +#include "grid_field.h" +#include "event.h" +#include "sphere.h" +#include "heap.h" + + +#define PI 3.141592653589793238462643 +#define SIZE 1.0 // size of box +#define VOLUMESPHERE pow(PI,((double)(DIM))/2.)/exp(lgamma(1+((double)(DIM))/2.)) // volume prefactor for sphere +#define DBL_EPSILON 2.2204460492503131e-016 // smallest # such that 1.0+DBL_EPSILON!=1.0 +#define M 1.0 + +//--------------------------------------------------------------------------- +// Class neighbor +//--------------------------------------------------------------------------- +class neighbor +{ + public: + int i; + + neighbor(int i_i); + + public: + virtual void Operation(int j, vector& pboffset) = 0; +}; + + +class box { + + public: + + // constructor and destructor + box(int N_i, double r_i, double growthrate_i, double maxpf_i, + double bidispersityratio, double bidispersityfraction, + double massratio, int hardwallBC); + ~box(); + + // Creating configurations + int Optimalngrids(); + int Optimalngrids2(double maxr); + void CreateSpheres(double temp); + void CreateSphere(int Ncurrent); + double Velocity(double temp); + void VelocityGiver(double temp); + void SetInitialEvents(); + void RecreateSpheres(const char* filename, double temp); + void ReadPositions(const char* filename); + void AssignCells(); + + // Predicting next event + event FindNextEvent(int i); + void CollisionChecker(event c); + event FindNextTransfer(int i); + event FindNextCollision(int i); + void ForAllNeighbors(int, vector, vector, neighbor&); + void PredictCollision(int i, int j, vector pboffset, + double& ctime, int& cpartner, + vector& cpartnerpboffset); + double CalculateCollision(int i, int j, vector pboffset); + double QuadraticFormula(double a, double b, double c); + + // Processing an event + void Process(int n); + void ProcessEvent(); + void Collision(event e); + void Transfer(event e); + void UpdateCell(int i, vector& celli); + void Synchronize(bool rescale); + void ChangeNgrids(int newngrids); + + // Debugging + void TrackPositions(); + void OutputEvents(); + void OutputCells(); + void GetInfo(); + int CheckSphereDiameters(); + + // Statistics + double Energy(); + double PackingFraction(); + void PrintStatistics(); + void RunTime(); + void WriteConfiguration(double* data); + + + //variables + + const int N; // number of spheres + + int ngrids; // number of cells in one direction + double maxpf; + double growthrate; // growth rate of the spheres + double r; // radius, defined at gtime = 0 + double gtime; // this is global clock + double rtime; // reset time, total time = rtime + gtime + double collisionrate; // average rate of collision between spheres + double bidispersityratio; // ratio of sphere radii + double bidispersityfraction; // fraction of smaller spheres + double massratio; // ratio of sphere masses + int hardwallBC; // =0 for periodic BC, =1 for hard wall + + + // statistics + double pressure; // pressure + double xmomentum; // exchanged momentum + double pf; // packing fraction + double energy; // kinetic energy + double energychange; + int ncollisions; // number of collisions + int ntransfers; // number of transfers + int nchecks; // number of checks + int neventstot; // total number of events + int ncycles; // counts # cycles for output + + time_t start, error, end; // run time of program + + // arrays + sphere *s; // array of spheres + grid_field cells; // array that keeps track of spheres in each cell + int *binlist; // linked-list for cells array + heap h; // event heap + vector *x; // positions of spheres.used for graphics +}; + + +//--------------------------------------------------------------------------- +// Predicts collisions, inherits neighbor operation +//--------------------------------------------------------------------------- +class collision : public neighbor +{ + public: + + box *b; + double ctime; + int cpartner; + vector cpartnerpboffset; + + public: + collision(int i_i, box *b); + + virtual void Operation(int j, vector& pboffset); +}; + +#endif diff --git a/src/spheres_poly/event.cpp b/src/spheres_poly/event.cpp new file mode 100644 index 0000000..c8c104d --- /dev/null +++ b/src/spheres_poly/event.cpp @@ -0,0 +1,65 @@ +#include "event.h" + +//============================================================== +//============================================================== +// Class Event +//============================================================== +//============================================================== + + +//============================================================== +// Constructor +//============================================================== +event::event(double time_i, int i_i, int j_i, vector v_i): + time(time_i), + i(i_i), + j(j_i), + v(v_i) +{ +} + +event::event(double time_i, int i_i, int j_i): + time(time_i), + i(i_i), + j(j_i) +{ +} + +event::event(const event& e) +{ + time = e.time; + i = e.i; + j = e.j; + v = e.v; +} + +event::event() +{ +} + + +//============================================================== +// Destructor +//============================================================== +event::~event() +{ +} + +void event::erase() +{ + time = dblINF; + i = 0; + j = 0; +} + +bool event::operator<(const event& e) const +{ + return e.time < time; +} + +bool event::operator>(const event& e) const +{ + return e.time > time; +} + + diff --git a/src/spheres_poly/event.h b/src/spheres_poly/event.h new file mode 100644 index 0000000..55dd1fc --- /dev/null +++ b/src/spheres_poly/event.h @@ -0,0 +1,58 @@ +#ifndef EVENT_H +#define EVENT_H + +#include "vector.h" +#define INF 100000000 +#define dblINF 100000000. + +class event { + + public: + + // constructor and destructor + event(double time_i, int i_i, int j_i, vector v_i); + event(double time_i, int i_i, int j_i); + event(const event& e); + event(); + + ~event(); + + bool operator<(const event&) const; + bool operator>(const event&) const; + void erase(); + + //variables + + + double time; // time of next collision + int i; // collision partner with lower number + int j; // collision partner with higher number + vector v; // virtual image + + /* 0<=j<=N binary collision between i and j + j=N+DIM+1+x transfer where x=-(k+1) for left wall + and x=k+1 for right wall + j=INF both check after event that did not altered motion of i and check after event that altered motion of i, i.e rescaling of velocities. I currently don't see need t o separate the two + + j=-1 check after collision + + Virtual identifiers as scalars...I think bad idea, but here's my work + there will be easily fixed problems if DIM>=10 + -x<=v<=x x=k+1, image in k direction + v=xy x,y + =-xy -x,y + =-yx x,-y + =yx -x,-y + v=xyz x,y,z + =-xyz -x,y,z + =-yxz x,-y,z + =-zxy x,y,-z + =zyx -x,-y,z + =xzy x,-y,-z + =yzx -x,y,-z + =-zyx -x,-y,-z + */ + +}; + +#endif diff --git a/src/spheres_poly/grid_field.h b/src/spheres_poly/grid_field.h new file mode 100644 index 0000000..10430ed --- /dev/null +++ b/src/spheres_poly/grid_field.h @@ -0,0 +1,148 @@ +/* + Packing of hard spheres via molecular dynamics + Developed by Monica Skoge, 2006, Princeton University + Contact: Aleksandar Donev (adonev@math.princeton.edu) with questions + This code may be used, modified and distributed freely. + Please cite: + + "Packing Hyperspheres in High-Dimensional Euclidean Spaces" + M. Skoge, A. Donev, F. H. Stillinger and S. Torquato, 2006 + + if you use these codes. +*/ + + +#ifndef GRID_FIELD_H +#define GRID_FIELD_H + +#include "vector.h" + +// ====================================================================== +// grid_field +// ====================================================================== + +// A field of V-vectors on a D dimensional manifold + +template +class grid_field { + + public: + int elements; + + private: + T* f; + vector size; // number of grid points for each dimension + vector offset; + + public: + + grid_field(); + grid_field(const vector&); + ~grid_field(); + + T& get(const vector&); + + vector get_size() const; + void set_size(const vector&); + void set_size(const int); + void initialize(const int i); +}; + + +// grid_field +// ~~~~~~~~~~~~ +template +grid_field::grid_field() + : f(0), elements(0) +{ +} + + +// grid_field +// ~~~~~~~~~~~~ +template +grid_field::grid_field(const vector& s) + : f(0) +{ + set_size(s); +} + + +// ~grid_field +// ~~~~~~~~~~~~~ +template +grid_field::~grid_field() +{ + if(f != 0) + delete[] f; +} + + +// get_size +// ~~~~~~~~ +template +inline vector grid_field::get_size() const +{ + return size; +} + + +// set_size +// ~~~~~~~~ +template +void grid_field::set_size(const vector& s) +{ + if(f != 0) + delete[] f; + + size = s; + + elements = 1; + for(int i=0; i +void grid_field::set_size(const int s) +{ + vector square; + + for(int k=0; k +inline T& grid_field::get(const vector& pos) +{ + int p=0; + for(int i=0; i +void grid_field::initialize(const int value) +{ + for(int i=0; i +#include "box.h" + + +//============================================================== +//============================================================== +// Class heap: Event heap used in hard spheres calculation +//============================================================== +//============================================================== + + +//============================================================== +// Constructor +//============================================================== +heap::heap(int maxsize_i): + maxsize(maxsize_i) +{ + a = new int[maxsize]; + index = new int[maxsize]; + + N = 0; // current number of events in heap +} + + +//============================================================== +// Constructor +//============================================================== +heap::heap(const heap &h) +{ + maxsize = h.maxsize; + a = h.a; + index = h.index; + N = h.N; // current number of events in heap + s = h.s; +} + + +//============================================================== +// Destructor +//============================================================== +heap::~heap() +{ + delete[] a; + delete[] index; +} + +//============================================================== +// Upheap +//============================================================== +void heap::upheap(int k) +{ + int i = a[k]; + + while ((k>1)&&(s[a[k/2]].nextevent.time > s[i].nextevent.time)) + { + a[k] = a[k/2]; + index[a[k/2]] = k; + k = k/2; + } + a[k] = i; + index[i] = k; +} + +//============================================================== +// Downheap +//============================================================== +void heap::downheap(int k) +{ + int j; + int i = a[k]; + + while(k <= N/2) + { + j = k+k; + if ((j < N)&&(s[a[j]].nextevent.time > s[a[j+1]].nextevent.time)) + j++; + if (s[i].nextevent.time <= s[a[j]].nextevent.time) + break; + a[k] = a[j]; + index[a[j]] = k; + k = j; + } + a[k] = i; + index[i] = k; +} + +//============================================================== +// Insert +//============================================================== +void heap::insert(int i) +{ + if (N >= maxsize) + std::cout << "error, N >= maxsize, cannot insert another event" << std::endl; + else + { + N++; + a[N] = i; + index[i] = N; + upheap(N); + } +} + + +//============================================================== +// Extract max +//============================================================== +int heap::extractmax() +{ + return a[1]; +} + +/* +//============================================================== +// Replace +//============================================================== +void heap::replace(int i) +{ + a[1] = i; + s[i].nextevent = t; + + if (!(e.time > s[a[1]].nextevent.time)) + { + if (!(s[a[1]].nextevent.j == INF))// must be check i'm changing to coll. at same time + std::cout << "error replaced non check with earlier time" << std::endl; + a[1] = i; + index[i] = 1; + } + else + { + a[0] = i; + downheap(0); + } +} +*/ + +//============================================================== +// Search +//============================================================== +int heap::search(int j) +{ + for (int k=1; k<=N; k++) + { + if (a[k] == j) + return k; + } + return -1; +} + +/* +//============================================================== +// Change +//============================================================== +void heap::change(int i) +{ + int iindex = index[i]; + + if (s[i].nextevent.time == s[a[iindex]].nextevent.time) + std::cout << "error changing an event to an equal time" << std::endl; + else if (s[i].nextevent.time > s[a[iindex]].nextevent.time) + std::cout << "error changing an event to a greater time" << std::endl; + + a[iindex] = i; + upheap(iindex); +} +*/ + +//============================================================== +// Print +//============================================================== +void heap::print() +{ + for (int k=1; k<=N; k++) + std::cout << k << " " << a[k] << " " << s[a[k]].nextevent.j << " " << s[a[k]].nextevent.time << std::endl; +} + + +//============================================================== +// Check index array +//============================================================== +void heap::checkindex() +{ + for (int k=1; k<=N; k++) + std::cout << k << " " << a[k] << " " << index[a[k]] << std::endl; +} diff --git a/src/spheres_poly/heap.h b/src/spheres_poly/heap.h new file mode 100644 index 0000000..91180e3 --- /dev/null +++ b/src/spheres_poly/heap.h @@ -0,0 +1,42 @@ +//--------------------------------------------------------------------------- +// Event heap maker +//--------------------------------------------------------------------------- + +#ifndef HEAP_H +#define HEAP_H + +#include "event.h" +#include "sphere.h" + +class heap { + + public: + + // constructor and destructor + heap(int maxsize); + heap(const heap &h); + ~heap(); + + // variables + int maxsize; // max allowed number of events + int N; // current number of events + int *a; + sphere *s; + int *index; // array of indices for each sphere + //event minevent; + + + // functions which operate on a binary heap + + void upheap(int k); + void downheap(int k); + void insert(int i); + void replace(int i); + int search(int j); + void change(int i); + int extractmax(); + void print(); + void checkindex(); + +}; +#endif diff --git a/src/spheres_poly/neighbor.cpp b/src/spheres_poly/neighbor.cpp new file mode 100644 index 0000000..920c099 --- /dev/null +++ b/src/spheres_poly/neighbor.cpp @@ -0,0 +1,40 @@ +#include "vector.h" +#include +#include "box.h" + + +//============================================================== +//============================================================== +// Class neighbor +//============================================================== +//============================================================== + +neighbor::neighbor(int i_i) + : i(i_i) +{ +} + + +//============================================================== +//============================================================== +// Class collision +//============================================================== +//============================================================== + +collision::collision(int i_i, box *b_i) + : neighbor(i_i), + b(b_i) +{ + ctime = dblINF; + cpartner = i; +} + + +//============================================================== +// Operation is finding the next collision from a given cell +//============================================================== +void collision::Operation(int j, vector& pboffset) +{ + b->PredictCollision(i, j, pboffset, ctime, cpartner, cpartnerpboffset); +} + diff --git a/src/spheres_poly/sphere.cpp b/src/spheres_poly/sphere.cpp new file mode 100644 index 0000000..0fbbd68 --- /dev/null +++ b/src/spheres_poly/sphere.cpp @@ -0,0 +1,66 @@ +#include "box.h" +#include +#include +#include +#include + +#include "vector.h" + +//============================================================== +//============================================================== +// Class Sphere: +//============================================================== +//============================================================== + + +//============================================================== +// Constructor +//============================================================== +sphere::sphere() +{ +} + + +//============================================================== +// Constructor +//============================================================== +sphere::sphere(const sphere& s) +{ + i = s.i; + x = s.x; + v = s.v; + cell = s.cell; + lutime = s.lutime; + nextevent = s.nextevent; + nextcollision = s.nextcollision; + r = s.r; + gr = s.gr; + m = s.m; + species=s.species; +} + +//============================================================== +// Constructor +//============================================================== +sphere::sphere(int i_i, vector x_i, vector cell_i, + double lutime_i, double r_i, double gr_i, double m_i, int species_i): + i(i_i), + x(x_i), + cell(cell_i), + lutime(lutime_i), + r(r_i), + gr(gr_i), + m(m_i), + species(species_i) +{ +} + +//============================================================== +// Destructor +//============================================================== +sphere::~sphere() +{ + +} + + diff --git a/src/spheres_poly/sphere.h b/src/spheres_poly/sphere.h new file mode 100644 index 0000000..28c64b4 --- /dev/null +++ b/src/spheres_poly/sphere.h @@ -0,0 +1,44 @@ +#ifndef SPHERE_H +#define SPHERE_H + + +#include "vector.h" + + +class sphere { + + public: + + // constructor and destructor + + sphere(); + sphere(const sphere& s); + sphere(int i_i, vector x, vector cell_i, double lutime_i, + double r_i, double gr_i, double m_i, int species_i); + ~sphere(); + + //variables + + int i; // sphere ID + + // impending event + event nextevent; // next event...can be collision or transfer + event nextcollision; // next collision if next event is transfer + // maybe nextnext event + + // past information + double lutime; // last update time + vector cell; // cell that it belongs to + vector x; // position + vector v; // velocity + double r; // sphere radius + double gr; // sphere growth rate + double m; // sphere mass + int species; // species number (not used during the MD) + // make sure efficent in memory + + + +}; + +#endif diff --git a/src/spheres_poly/spheres.cpp b/src/spheres_poly/spheres.cpp new file mode 100644 index 0000000..9b14b01 --- /dev/null +++ b/src/spheres_poly/spheres.cpp @@ -0,0 +1,64 @@ +//=========================================================== +//=========================================================== +//=========================================================== +// +// Molecular dynamics simulation of hardspheres +// +//=========================================================== +//=========================================================== +//=========================================================== + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "box.h" +#include "sphere.h" +#include "event.h" +#include "heap.h" + + +double *spherespp(int N, double bidispersityratio, double bidispersityfraction, double maxpf, double maxpressure, double maxcollisionrate) +{ + int eventspercycle = 20; // # events per sphere per cycle + double initialpf = 0.01; // initial packing fraction + double temp = 0.2; // initial temp., use 0 for zero velocities + double growthrate = 0.001; // growth rate + double massratio = 1.; // ratio of sphere masses + int hardwallBC = 0; // 0 for periodic, 1 for hard wall BC + + double d, r; // initial diameter and radius of spheres + + r = pow(initialpf*pow(SIZE, DIM)/(N*VOLUMESPHERE), 1.0/((double)(DIM))); + + box b(N, r, growthrate, maxpf, bidispersityratio, + bidispersityfraction, massratio, hardwallBC); + + b.CreateSpheres(temp); + + + while ((b.collisionrate < maxcollisionrate) && (b.pf < maxpf) && (b.pressure < maxpressure)) + { + b.Process(eventspercycle*N); + b.Synchronize(true); + } + + double *data = (double *)malloc(DIM * N * sizeof(double)); + + b.WriteConfiguration(data); + + return data; +} + +extern "C" { + double *spheres(int N, double bidispersityratio, double bidispersityfraction, double maxpf, double maxpressure, double maxcollisionrate) { + return spherespp(N, bidispersityratio, bidispersityfraction, maxpf, maxpressure, maxcollisionrate); + } +} + + diff --git a/src/spheres_poly/spheres.h b/src/spheres_poly/spheres.h new file mode 100644 index 0000000..44d8052 --- /dev/null +++ b/src/spheres_poly/spheres.h @@ -0,0 +1,4 @@ + +#pragma once + +double *spheres(int N, double bidispersityratio, double bidispersityfraction, double maxpf, double maxpressure, double maxcollisionrate); diff --git a/src/spheres_poly/vector.h b/src/spheres_poly/vector.h new file mode 100644 index 0000000..9ee481f --- /dev/null +++ b/src/spheres_poly/vector.h @@ -0,0 +1,471 @@ +#ifndef VECTOR_H +#define VECTOR_H + +#include +#include + +#define DIM 2 + +template +class vector { + + public: + T x[D]; + + public: + + vector(); + vector(const T[D]); + vector(const vector&); + ~vector(); + + vector& operator+=(const vector&); + vector& operator-=(const vector&); + vector& operator*=(const T); + vector& operator/=(const T); + vector operator+(const vector&) const; + vector operator-(const vector&) const; + vector operator*(const T) const; + vector operator/(const T) const; + vector operator%(const T) const; + bool operator==(const vector &a) const; + + vector integer() const; + vector Double() const; + static vector integer(const vector&); + static vector Double(const vector&); + + T& operator[](const unsigned int); + + double dot(const vector&) const; + static double dot(const vector&, const vector&); + + double norm_squared() const; + static double norm_squared(const vector&); + + void read(std::ifstream&); + void write(std::ofstream&) const; +}; + + +template +std::ostream& operator<<(std::ostream&, const vector&); + + +// constructor +// ~~~~~~~~~~~ +template +vector::vector() +{ + for(int k=0; k +vector::vector(const T x_i[D]) +{ + for(int k=0; k +vector::vector(const vector &v) +{ + for(int k=0; k +vector::~vector() +{ +} + + +// += +// ~~ +template +inline vector& vector::operator+=(const vector &v) +{ + for(int k=0; k +inline vector& vector::operator-=(const vector &v) +{ + for(int k=0; k +inline vector& vector::operator*=(const T s) +{ + for(int k=0; k +inline vector& vector::operator/=(const T s) +{ + for(int k=0; k +inline vector vector::operator+(const vector &a) const +{ + vector c; + + for(int k=0; k +inline vector vector::operator-(const vector &a) const +{ + vector c; + + for(int k=0; k +inline vector vector::operator*(const T s) const +{ + vector c; + + for(int k=0; k +inline vector vector::operator/(const T s) const +{ + vector c; + + for(int k=0; k +inline bool vector::operator==(const vector &a) const +{ + for(int k=0; k +inline vector vector::operator%(const T s) const +{ + vector c; + + for(int k=0; k +inline vector vector::integer() const +{ + vector c; + + for(int k=0; k +inline vector vector::integer(const vector& v) +{ + return v.integer(); +} + + +// double +// ~~~~~~~ +template +inline vector vector::Double() const +{ + vector c; + + for(int k=0; k +inline vector vector::Double(const vector& v) +{ + return v.Double(); +} + + + +// [] +// ~~ +template +inline T& vector::operator[](const unsigned int i) +{ + return x[i]; +} + + +// Dot +// ~~~ +template +inline double vector::dot(const vector &a) const +{ + double d=0; + + for(int k=0; k +inline double vector::dot(const vector &a, const vector &b) +{ + return a.dot(b); +} + + +// NormSquared +// ~~~~~~~~~~~ +template +inline double vector::norm_squared() const +{ + return dot(*this, *this); +} + +template +inline double vector::norm_squared(const vector& v) +{ + return v.norm_squared(); +} + + +// read +// ~~~~ +template +void vector::read(std::ifstream& in) +{ + in.read((char*)x, sizeof(T)*D); +} + +// write +// ~~~~~ +template +void vector::write(std::ofstream& out) const +{ + out.write((const char*)x, sizeof(T)*D); +} + + + +// Insertion +// ~~~~~~~~~ +template +std::ostream& operator<<(std::ostream& os, const vector& v) +{ + os << "("; + + for(int k=0; k +class vector_field { + + public: + int elements; + + private: + vector* f; + vector size; // number of grid points for each dimension + vector offset; + + public: + + vector_field(); + vector_field(const vector&); + ~vector_field(); + + vector get_size() const; + void set_size(const vector&); + + vector& get(const vector&); + + void read(std::ifstream&); + void write(std::ofstream&) const; + + static void swap(vector_field&, vector_field&); +}; + + +// vector_field +// ~~~~~~~~~~~~ +template +vector_field::vector_field() + : f(0), elements(0) +{ +} + + +// vector_field +// ~~~~~~~~~~~~ +template +vector_field::vector_field(const vector& s) + : f(0) +{ + set_size(s); +} + +// ~vector_field +// ~~~~~~~~~~~~~ +template +vector_field::~vector_field() +{ + if(f != 0) + delete[] f; +} + +// get_size +// ~~~~~~~~ +template +inline vector vector_field::get_size() const +{ + return size; +} + +// set_size +// ~~~~~~~~ +template +void vector_field::set_size(const vector& s) +{ + if(f != 0) + delete[] f; + + size = s; + + elements = 1; + for(int i=0; i[elements]; +} + +// get +// ~~~ +template +inline vector& vector_field::get(const vector& pos) +{ + int p=0; + for(int i=0; i +void vector_field::read(std::ifstream& in) +{ + in.read((char*)f, elements*sizeof(T)*V); +} + + +// write +// ~~~~~ +template +void vector_field::write(std::ofstream& out) const +{ + out.write((const char*)f, elements*sizeof(T)*V); +} + +// swap +// ~~~~ +template +void vector_field::swap(vector_field& v1, + vector_field& v2) +{ + vector* f; + + f = v1.f; + v1.f = v2.f; + v2.f = f; +} + + + +#endif -- cgit v1.2.3-70-g09d2