summaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/bound_set.c102
-rw-r--r--src/break_edge.c51
-rw-r--r--src/corr_test.c4
-rw-r--r--src/correlations.c299
-rw-r--r--src/data.c35
-rw-r--r--src/factor_update.c68
-rw-r--r--src/fitting_functions.c290
-rw-r--r--src/fortune/defs.h134
-rw-r--r--src/fortune/edgelist.c136
-rw-r--r--src/fortune/geometry.c184
-rw-r--r--src/fortune/heap.c94
-rw-r--r--src/fortune/main.c287
-rw-r--r--src/fortune/memory.c44
-rw-r--r--src/fortune/output.c46
-rw-r--r--src/fortune/voronoi.c103
-rw-r--r--src/fracture.c6
-rw-r--r--src/fracture.h91
-rw-r--r--src/frame_create.c167
-rw-r--r--src/gen_laplacian.c137
-rw-r--r--src/gen_voltcurmat.c36
-rw-r--r--src/geometry.c55
-rw-r--r--src/get_dual_clusters.c36
-rw-r--r--src/graph_components.c185
-rw-r--r--src/graph_create.c299
-rw-r--r--src/graph_free.c19
-rw-r--r--src/graph_genfunc.c19
-rw-r--r--src/net.c141
-rw-r--r--src/net_conductivity.c35
-rw-r--r--src/net_currents.c51
-rw-r--r--src/net_fracture.c67
-rw-r--r--src/net_voltages.c40
-rw-r--r--src/rand.c23
-rw-r--r--src/spheres_poly/box.cpp1157
-rw-r--r--src/spheres_poly/box.h169
-rw-r--r--src/spheres_poly/event.cpp65
-rw-r--r--src/spheres_poly/event.h58
-rw-r--r--src/spheres_poly/grid_field.h148
-rw-r--r--src/spheres_poly/heap.cpp186
-rw-r--r--src/spheres_poly/heap.h42
-rw-r--r--src/spheres_poly/neighbor.cpp40
-rw-r--r--src/spheres_poly/sphere.cpp66
-rw-r--r--src/spheres_poly/sphere.h44
-rw-r--r--src/spheres_poly/spheres.cpp64
-rw-r--r--src/spheres_poly/spheres.h4
-rw-r--r--src/spheres_poly/vector.h471
45 files changed, 13 insertions, 5785 deletions
diff --git a/src/bound_set.c b/src/bound_set.c
deleted file mode 100644
index 65f1e6f..0000000
--- a/src/bound_set.c
+++ /dev/null
@@ -1,102 +0,0 @@
-
-#include "fracture.h"
-
-double th_p(double x, double y, double th) {
- if (x >= 0 && y >= 0) return th;
- else if (x < 0 && y >= 0) return M_PI - th;
- else if (x < 0 && y < 0) return th - M_PI;
- else return -th;
-
-}
-
-double u_y(double x, double y) {
- double r = sqrt(pow(x, 2) + pow(y, 2));
- double th = th_p(x, y, atan(fabs(y / x)));
-
- return sqrt(r) * sin(th / 2);
-}
-
-void bound_set_embedded(double *bound, const graph_t *g, double notch_len) {
- uint_t L = g->L;
-
- for (uint_t i = 0; i < L / 2; i++) {
- double x1, y1, x2, y2, x3, y3, x4, y4;
- x1 = (2. * i + 1.) / L - notch_len; y1 = 0.5 - 1.;
- x2 = (2. * i + 1.) / L - notch_len; y2 = 0.5 - 0.;
- y3 = (2. * i + 1.) / L - 0.5; x3 = 0.5 - 1.;
- y4 = (2. * i + 1.) / L - 0.5; x4 = 0.5 - 0.;
-
- bound[g->b[g->bi[0] + i]] = u_y(x1, y1);
- bound[g->b[g->bi[1] + i]] = u_y(x2, y2);
- bound[g->b[g->bi[2] + i]] = u_y(x3, y3);
- bound[g->b[g->bi[3] + i]] = u_y(x4, y4);
- }
-}
-
-bool is_in(uint_t len, uint_t *list, uint_t element) {
- for (uint_t i = 0; i < len; i++) {
- if (list[i] == element) {
- return true;
- }
- }
- return false;
-}
-
-cholmod_dense *bound_set(const graph_t *g, bool vb, double notch_len, cholmod_common *c) {
-
- uint_t dim = g->nv;
-
- if (vb && g->boundary != TORUS_BOUND) {
- dim -= g->bi[g->nb];
- } else if (!vb) {
- dim += 2;
- }
-
- cholmod_dense *boundary = CHOL_F(zeros)(dim, 1, CHOLMOD_REAL, c);
- double *bound = (double *)boundary->x;
-
- switch (g->boundary) {
- case TORUS_BOUND:
- for (uint_t i = 0; i < g->bi[1]; i++) {
- uint_t be = g->b[i];
- uint_t v1 = g->ev[2 * be];
- uint_t v2 = g->ev[2 * be + 1];
- double v1y = g->vx[2 * v1 + 1];
- double v2y = g->vx[2 * v2 + 1];
-
- uint_t ind = v1y < v2y ? 0 : 1;
-
- bound[g->ev[2 * be + ind]] += 1;
- bound[g->ev[2 * be + !ind]] -= 1;
- }
- break;
- /*
- case EMBEDDED_BOUND:
- bound_set_embedded(bound, g, notch_len);
- break;
- */
- default:
- if (vb) {
- for (uint_t i = 0; i < dim; i++) {
- uint_t v = g->nbi[i];
- for (uint_t j = 0; j < g->vei[v+1] - g->vei[v]; j++) {
- uint_t e = g->ve[g->vei[v] + j];
- uint_t v0 = g->ev[2 * e];
- uint_t v1 = g->ev[2 * e + 1];
-
- if (g->bq[v0] || g->bq[v1]) {
- uint_t vv = v0 == v ? v1 : v0;
- if (is_in(g->bi[1], g->b, vv)) {
- bound[i]++;
- }
- }
- }
- }
- } else {
- bound[g->nv] = 1;
- bound[g->nv + 1] = -1;
- }
- }
-
- return boundary;
-}
diff --git a/src/break_edge.c b/src/break_edge.c
deleted file mode 100644
index 2f112c2..0000000
--- a/src/break_edge.c
+++ /dev/null
@@ -1,51 +0,0 @@
-
-#include "fracture.h"
-
-
-bool break_edge(net_t *net, uint_t edge, cholmod_common *c) {
- net->fuses[edge] = true;
- net->num_broken++;
-
- double *x = (double *)net->boundary_cond->x;
- const graph_t *g = net->graph;
-
- uint_t v1 = net->graph->ev[2 * edge];
- uint_t v2 = net->graph->ev[2 * edge + 1];
-
- uint_t s1 = v1 < v2 ? v1 : v2;
- uint_t s2 = v1 < v2 ? v2 : v1;
-
- if (net->graph->boundary == TORUS_BOUND) {
- if (net->graph->bq[edge]) {
- double v1y = g->vx[2 * v1 + 1];
- double v2y = g->vx[2 * v2 + 1];
- uint_t ind = v1y < v2y ? 0 : 1;
- x[g->ev[2 * edge + ind]] -= 1;
- x[g->ev[2 * edge + !ind]] += 1;
- }
- if (net->factor != NULL) {
- factor_update(net->factor, s1, s2, c);
- }
- } else if (net->voltage_bound) {
- if (g->bq[v1] || g->bq[v2]) {
- uint_t vv = g->bq[v1] ? v2 : v1;
- uint_t uu = v1 == vv ? v2 : v1;
- if (is_in(g->bi[1], g->b, uu)) {
- x[g->bni[vv]] -= 1;
- }
- if (net->factor != NULL) {
- factor_update2(net->factor, g->bni[vv], c);
- }
- } else {
- if (net->factor != NULL) {
- factor_update(net->factor, g->bni[s1], g->bni[s2], c);
- }
- }
- } else {
- if (net->factor != NULL) {
- factor_update(net->factor, s1, s2, c);
- }
- }
-
- return true;
-}
diff --git a/src/corr_test.c b/src/corr_test.c
index b6f8a18..01b3e11 100644
--- a/src/corr_test.c
+++ b/src/corr_test.c
@@ -7,7 +7,7 @@ int main() {
unsigned int width = 64;
- graph_t *network = graph_create(VORONOI_LATTICE, TORUS_BOUND, 128, false, &c);
+ graph_t *network = graph_create(VORONOI_LATTICE, TORUS_BOUND, 128, false);
net_t *instance = net_create(network, 1e-14, 3, 0, true, &c);
net_fracture(instance, &c, 1e-10);
@@ -19,7 +19,7 @@ int main() {
printf("\n");
net_free(instance, &c);
- graph_free(network, &c);
+ graph_free(network);
CHOL_F(finish)(&c);
diff --git a/src/correlations.c b/src/correlations.c
deleted file mode 100644
index f403ab8..0000000
--- a/src/correlations.c
+++ /dev/null
@@ -1,299 +0,0 @@
-
-#include "fracture.h"
-
-// I implemented a fibonacci heap because wikipedia said it would be fast
-
-struct fibn {
- uint_t value;
- uint_t priority;
- uint_t rank;
- bool marked;
- struct fibn *first_child;
- struct fibn *parent;
- struct fibn *prev;
- struct fibn *next;
-};
-
-typedef struct {
- struct fibn *min;
- uint_t rank;
- uint_t trees;
- struct fibn **vertex_to_node;
-} fibh;
-
-uint_t fib_findmin(fibh *heap) {
- return heap->min->value;
-}
-
-void ll_setparent(struct fibn *ll, struct fibn *parent) {
- struct fibn *curnode = ll->next;
- ll->parent = parent;
- while (curnode != ll) {
- curnode->parent = parent;
- curnode = curnode->next;
- }
-}
-
-struct fibn *ll_merge(struct fibn *ll1, struct fibn *ll2) {
- if (ll1 == NULL)
- return ll2;
- if (ll2 == NULL)
- return ll1;
-
- // link the beginning of list one to the end of list two and vice versa
- struct fibn *ll1_beg = ll1;
- struct fibn *ll1_end = ll1->prev;
- struct fibn *ll2_beg = ll2;
- struct fibn *ll2_end = ll2->prev;
-
- ll1_beg->prev = ll2_end;
- ll1_end->next = ll2_beg;
- ll2_beg->prev = ll1_end;
- ll2_end->next = ll1_beg;
-
- return ll1;
-}
-
-struct fibn *ll_delroot(struct fibn *ll) {
- if (ll == NULL)
- return NULL;
-
- if (ll->next == ll) {
- return NULL;
- }
-
- struct fibn *ll_beg = ll->next;
- struct fibn *ll_end = ll->prev;
-
- ll_beg->prev = ll_end;
- ll_end->next = ll_beg;
-
- ll->next = ll;
- ll->prev = ll;
-
- return ll_beg;
-}
-
-void fib_insert(fibh *heap, uint_t value, uint_t priority) {
- struct fibn *newnode = calloc(1, sizeof(struct fibn));
- newnode->value = value;
- newnode->priority = priority;
- newnode->next = newnode;
- newnode->prev = newnode;
-
- heap->min = ll_merge(heap->min, newnode);
-
- if (priority < heap->min->priority) {
- heap->min = newnode;
- }
-
- heap->vertex_to_node[value] = newnode;
-
- heap->trees++;
-}
-
-void fib_deletemin(fibh *heap) {
- uint_t min_rank = heap->min->rank;
- struct fibn *min_children = heap->min->first_child;
-
- struct fibn *trees = ll_delroot(heap->min);
- heap->vertex_to_node[heap->min->value] = NULL;
- free(heap->min);
-
- if (min_children != NULL)
- ll_setparent(min_children, NULL);
- trees = ll_merge(trees, min_children);
-
- heap->trees += min_rank - 1;
-
- if (trees == NULL) {
- // min had no children and was only tree, return empty heap
- heap->min = NULL;
- heap->rank = 0;
- } else {
- // min had children or there were other trees
- uint_t min_val = UINT_MAX;
- struct fibn *min_ptr = NULL;
- uint_t max_rank = 0;
- struct fibn *curnode = trees;
- for (uint_t i = 0; i < heap->trees; i++) {
- if (curnode->priority < min_val) {
- min_ptr = curnode;
- min_val = curnode->priority;
- }
- if (curnode->rank > max_rank)
- max_rank = curnode->rank;
- curnode = curnode->next;
- }
-
- if (min_ptr == NULL)
- min_ptr = trees;
- heap->min = min_ptr;
- heap->rank = max_rank;
-
- struct fibn **rankslots = calloc(max_rank + 1, sizeof(struct fibn *));
- curnode = heap->min;
- while (curnode != rankslots[curnode->rank]) {
- if (rankslots[curnode->rank] == NULL) {
- rankslots[curnode->rank] = curnode;
- curnode = curnode->next;
- } else {
- struct fibn *oldnode = rankslots[curnode->rank];
- rankslots[curnode->rank] = NULL;
- struct fibn *smaller =
- curnode->priority < oldnode->priority || curnode == heap->min
- ? curnode
- : oldnode;
- struct fibn *larger = smaller == curnode ? oldnode : curnode;
- ll_delroot(larger);
- ll_setparent(larger, smaller);
- struct fibn *smaller_children = smaller->first_child;
- smaller->first_child = ll_merge(smaller_children, larger);
- heap->trees--;
- smaller->rank++;
- if (smaller->rank > heap->rank) {
- heap->rank = smaller->rank;
- rankslots =
- realloc(rankslots, (heap->rank + 1) * sizeof(struct fibn *));
- rankslots[heap->rank] = NULL;
- }
- curnode = smaller;
- }
- }
- free(rankslots);
- }
-}
-
-void fib_decreasekey(fibh *heap, uint_t value, uint_t new_priority) {
- struct fibn *node = heap->vertex_to_node[value];
- if (node != NULL) {
- node->priority = new_priority;
- if (node->parent != NULL) {
- if (node->priority < node->parent->priority) {
-
- struct fibn *curnode = node;
- curnode->marked = true;
- while (curnode->parent != NULL && curnode->marked) {
- struct fibn *oldparent = curnode->parent;
- oldparent->rank--;
- oldparent->first_child = ll_delroot(curnode);
- ll_setparent(curnode, NULL);
- ll_merge(heap->min, curnode);
- heap->trees++;
- if (curnode->marked) {
- curnode->marked = false;
- }
-
- if (oldparent->marked) {
- curnode = oldparent;
- } else {
- oldparent->marked = true;
- break;
- }
- }
- }
- }
-
- if (new_priority < heap->min->priority) {
- heap->min = node;
- }
- }
-}
-
-uint_t *dijkstra(const graph_t *network, uint_t source) {
- uint_t nv = network->dnv;
- uint_t *vei = network->dvei;
- uint_t *ve = network->dve;
- uint_t *ev = network->dev;
-
- uint_t *dist = (uint_t *)calloc(nv, sizeof(uint_t));
- fibh *Q = (fibh *)calloc(1, sizeof(fibh));
- Q->vertex_to_node = (struct fibn **)calloc(nv, sizeof(struct fibn *));
-
- for (uint_t i = 0; i < nv; i++) {
- if (i != source) {
- dist[i] = UINT_MAX;
- }
-
- fib_insert(Q, i, dist[i]);
- }
-
- while (Q->min != NULL) {
- uint_t u = fib_findmin(Q);
- fib_deletemin(Q);
-
- for (uint_t i = 0; i < vei[u + 1] - vei[u]; i++) {
- uint_t e = ve[vei[u] + i];
- uint_t v = ev[2 * e] == u ? ev[2 * e + 1] : ev[2 * e];
- uint_t alt = dist[u] + 1;
- if (alt < dist[v]) {
- dist[v] = alt;
- fib_decreasekey(Q, v, alt);
- }
- }
- }
-
- free(Q->vertex_to_node);
- free(Q);
-
- return dist;
-}
-
-uint_t **get_dists(const graph_t *network) {
- uint_t nv = network->dnv;
- uint_t **dists = (uint_t **)malloc(nv * sizeof(uint_t *));
-
- for (uint_t i = 0; i < nv; i++) {
- dists[i] = dijkstra(network, i);
- }
-
- return dists;
-}
-
-double *get_corr(net_t *instance, uint_t **dists, cholmod_common *c) {
- uint_t nv = instance->graph->dnv;
- uint_t ne = instance->graph->ne;
- uint_t *ev = instance->graph->dev;
- bool nulldists = false;
- if (dists == NULL) {
- dists = get_dists(instance->graph);
- nulldists = true;
- }
- double *corr = calloc(nv, sizeof(double));
- uint_t *marks = get_clusters(instance, c);
- uint_t *numat = calloc(nv, sizeof(uint_t));
-
- for (uint_t i = 0; i < ne; i++) {
- uint_t v1 = ev[2 * i];
- uint_t v2 = ev[2 * i + 1];
- for (uint_t j = 0; j < ne; j++) {
- uint_t v3 = ev[2 * j];
- uint_t v4 = ev[2 * j + 1];
- uint_t dist1 = dists[v1][v3];
- uint_t dist2 = dists[v1][v4];
- uint_t dist3 = dists[v2][v3];
- uint_t dist4 = dists[v2][v4];
- uint_t dist = (dist1 + dist2 + dist3 + dist4) / 4;
- corr[dist] += instance->fuses[i] && instance->fuses[j];
- numat[dist]++;
- }
- }
-
- for (uint_t i = 0; i < nv; i++) {
- if (numat[i] > 0) {
- corr[i] /= numat[i];
- }
- }
-
- if (nulldists) {
- for (int i = 0; i < nv; i++) {
- free(dists[i]);
- }
- free(dists);
- }
-
- free(marks);
-
- return corr;
-}
-
diff --git a/src/data.c b/src/data.c
deleted file mode 100644
index 2047c44..0000000
--- a/src/data.c
+++ /dev/null
@@ -1,35 +0,0 @@
-
-#include "fracture.h"
-
-data_t *data_create(uint_t ne) {
- data_t *data = malloc(1 * sizeof(data_t));
- assert(data != NULL);
-
- data->num_broken = 0;
-
- data->break_list = (uint_t *)malloc(ne * sizeof(uint_t));
- assert(data->break_list != NULL);
-
- data->extern_field = (long double *)malloc(ne * sizeof(long double));
- assert(data->extern_field != NULL);
-
- data->conductivity = (double *)malloc(ne * sizeof(double));
- assert(data->conductivity != NULL);
-
- return data;
-}
-
-void data_free(data_t *data) {
- free(data->break_list);
- free(data->extern_field);
- free(data->conductivity);
- free(data);
-}
-
-void data_update(data_t *data, uint_t last_broke, long double strength, double conductivity) {
- data->break_list[data->num_broken] = last_broke;
- data->extern_field[data->num_broken] = strength;
- data->conductivity[data->num_broken] = conductivity;
- data->num_broken++;
-}
-
diff --git a/src/factor_update.c b/src/factor_update.c
deleted file mode 100644
index a51705a..0000000
--- a/src/factor_update.c
+++ /dev/null
@@ -1,68 +0,0 @@
-
-#include "fracture.h"
-
-void factor_update(cholmod_factor *factor, uint_t v1, uint_t v2, cholmod_common *c) {
- uint_t n = factor->n;
-
- cholmod_sparse *update_mat =
- CHOL_F(allocate_sparse)(n, n, 2, true, true, 0, CHOLMOD_REAL, c);
-
- uint_t s1, s2;
- s1 = v1 < v2 ? v1 : v2;
- s2 = v1 > v2 ? v1 : v2;
-
- int_t *pp = (int_t *)update_mat->p;
- int_t *ii = (int_t *)update_mat->i;
- double *xx = (double *)update_mat->x;
-
- for (uint_t i = 0; i <= s1; i++) {
- pp[i] = 0;
- }
-
- for (uint_t i = s1 + 1; i <= n; i++) {
- pp[i] = 2;
- }
-
- ii[0] = s1;
- ii[1] = s2;
- xx[0] = 1;
- xx[1] = -1;
-
- cholmod_sparse *perm_update_mat = CHOL_F(submatrix)(
- update_mat, factor->Perm, factor->n, NULL, -1, true, true, c);
-
- CHOL_F(updown)(false, perm_update_mat, factor, c);
-
- CHOL_F(free_sparse)(&perm_update_mat, c);
- CHOL_F(free_sparse)(&update_mat, c);
-}
-
-void factor_update2(cholmod_factor *factor, uint_t v, cholmod_common *c) {
- uint_t n = factor->n;
-
- cholmod_sparse *update_mat =
- CHOL_F(allocate_sparse)(n, n, 1, true, true, 0, CHOLMOD_REAL, c);
-
- int_t *pp = (int_t *)update_mat->p;
- int_t *ii = (int_t *)update_mat->i;
- double *xx = (double *)update_mat->x;
-
- for (uint_t i = 0; i <= v; i++) {
- pp[i] = 0;
- }
-
- for (uint_t i = v + 1; i <= n; i++) {
- pp[i] = 1;
- }
-
- ii[0] = v;
- xx[0] = 1;
-
- cholmod_sparse *perm_update_mat = CHOL_F(submatrix)(
- update_mat, factor->Perm, factor->n, NULL, -1, true, true, c);
-
- CHOL_F(updown)(false, perm_update_mat, factor, c);
-
- CHOL_F(free_sparse)(&perm_update_mat, c);
- CHOL_F(free_sparse)(&update_mat, c);
-}
diff --git a/src/fitting_functions.c b/src/fitting_functions.c
deleted file mode 100644
index aa7a305..0000000
--- a/src/fitting_functions.c
+++ /dev/null
@@ -1,290 +0,0 @@
-
-#include "fracture.h"
-#include <gsl/gsl_sf_erf.h>
-#include <gsl/gsl_sf_laguerre.h>
-#include <gsl/gsl_matrix.h>
-#include <gsl/gsl_vector.h>
-#include <gsl/gsl_blas.h>
-#include <gsl/gsl_multifit_nlinear.h>
-
-struct model_params {
- double nu;
- double ac;
- double mc;
- double ac02;
- double ac03;
- double ac12;
- double ac13;
- double bc02;
- double bc03;
- double bc12;
- double bc13;
- double Ac20;
- double Ac21;
- double Ac22;
- double Ac23;
- double Ac30;
- double Ac31;
- double Ac32;
- double Ac33;
- double Bc20;
- double Bc21;
- double Bc22;
- double Bc23;
- double Bc30;
- double Bc31;
- double Bc32;
- double Bc33;
-};
-
-struct data {
- uint_t nc;
- uint_t *Ls_c;
- double *betas_c;
- double *vals_c2;
- double *vals_c3;
-};
-
-double Jc(uint_t n, uint_t o, struct model_params *par, double x) {
- double nu, ac, mc, ac0n, bc0n, *Acn, yc, sum;
-
- nu = par->nu;
- ac = par->ac;
- mc = par->mc;
- ac0n = *(&(par->ac02) + (n - 2) * sizeof(double));
- bc0n = *(&(par->bc02) + (n - 2) * sizeof(double));
- Acn = &(par->Ac20) + (n - 2) * o * sizeof(double);
-
- yc = (log(x) - mc) / ac;
-
- sum = 0;
-
- for (uint_t i = 0; i < o; i++) {
- sum += Acn[i] * gsl_sf_laguerre_n(i, 0, yc);
- }
-
- return exp(ac0n + bc0n * (2 - gsl_sf_erf(yc)) + exp(-gsl_pow_2(yc)) * sum);
-}
-
-double Kc(uint_t n, uint_t o, struct model_params *par, double x) {
- double nu, ac, mc, ac1n, bc1n, *Bcn, yc, sum;
-
- nu = par->nu;
- ac = par->ac;
- mc = par->mc;
- ac1n = *(&(par->ac12) + (n - 2) * sizeof(double));
- bc1n = *(&(par->bc12) + (n - 2) * sizeof(double));
- Bcn = &(par->Bc20) + (n - 2) * o * sizeof(double);
-
- yc = (log(x) - mc) / ac;
-
- sum = 0;
-
- for (uint_t i = 0; i < o; i++) {
- sum += Bcn[i] * gsl_sf_laguerre_n(i, 0, yc);
- }
-
- return exp(ac1n + bc1n * (2 - gsl_sf_erf(yc)) + exp(-gsl_pow_2(yc)) * sum);
-}
-
-double sc(uint_t n, uint_t o, struct model_params *par, uint_t L, double beta) {
- double nu, bLnu, deltanu, sigmanu, tau, Lntau, Ldelta;
-
- nu = par->nu;
- deltanu = 1.5;
- sigmanu = 48. / 91;
- tau = 187. / 91;
-
- bLnu = beta * exp(log(L) / nu);
- Lntau = exp((n + 1 - tau) / sigmanu * log(L));
- Ldelta = exp(-deltanu * log(L));
-
- return Lntau * (Jc(n, o, par, bLnu) + Ldelta * Kc(n, 0, par, bLnu));
-}
-
-int f_moms(const gsl_vector *x, void *params, gsl_vector *f) {
- struct data *dat = (struct data *)params;
-
- for (uint_t i = 0; i < dat->nc; i++) {
- double f2i = sc(2, 4, (struct model_params *)x->data, dat->Ls_c[i], dat->betas_c[i]);
- double F2i = dat->vals_c2[i];
-
- f2i = log(f2i);
- F2i = log(F2i);
-
- gsl_vector_set(f, i, f2i - F2i);
-
- double f3i = sc(3, 4, (struct model_params *)x->data, dat->Ls_c[i], dat->betas_c[i]);
- double F3i = dat->vals_c3[i];
-
- f3i = log(f3i);
- F3i = log(F3i);
-
- gsl_vector_set(f, dat->nc + i, f3i - F3i);
- }
-
- return GSL_SUCCESS;
-}
-
-void mom(uint_t len, uint_t n, uint32_t *data, double result[2]) {
- uint_t total = 0;
- double mom = 0;
- double mom_err = 0;
-
- for (uint_t i = 0; i < len; i++) {
- uint32_t datai = data[i];
- double in = pow(i, n);
-
- total += datai;
- mom += in * datai;
- mom_err += gsl_pow_2(in) * datai;
- }
-
- double momf = mom / total;
- double momf_err = momf * sqrt(mom_err / gsl_pow_2(mom) + 1 / total);
-
- result[0] = momf;
- result[1] = momf_err;
-}
-
-void
-callback(const size_t iter, void *params,
- const gsl_multifit_nlinear_workspace *w)
-{
- gsl_vector *f = gsl_multifit_nlinear_residual(w);
-
- fprintf(stderr, "iter %2zu: |f(x)| = %.4f\n",
- iter,
- gsl_blas_dnrm2(f));
-}
-
-
-int main(int argc, char *argv[]) {
- struct data *d = (struct data *)malloc(sizeof(struct data));
- uint_t nc = argc - 1;
- uint_t *Ls_c = (uint_t *)malloc(nc * sizeof(uint_t));
- double *betas_c = (double *)malloc(nc * sizeof(double));
- double *vals_c2 = (double *)malloc(nc * sizeof(double));
- double *weights_c2 = (double *)malloc(nc * sizeof(double));
- double *vals_c3 = (double *)malloc(nc * sizeof(double));
- double *weights_c3 = (double *)malloc(nc * sizeof(double));
-
- double chisq, chisq0;
-
- d->nc = nc;
- d->Ls_c = Ls_c;
- d->betas_c = betas_c;
- d->vals_c2 = vals_c2;
- d->vals_c3 = vals_c3;
-
- for (uint_t i = 0; i < nc; i++) {
- char *fn = argv[1 + i];
- char *buff = malloc(20 * sizeof(char));
- uint_t pos = 0; uint_t c = 0;
- while (c < 4) {
- if (fn[pos] == '_') {
- c++;
- }
- pos++;
- }
- c = 0;
- while (fn[pos] != '_') {
- buff[c] = fn[pos];
- pos++;
- c++;
- }
- buff[c] = '\0';
- Ls_c[i] = atoi(buff);
- c = 0;
- pos++;
- while (fn[pos] != '_') {
- buff[c] = fn[pos];
- pos++;
- c++;
- }
- buff[c] = '\0';
- betas_c[i] = atof(buff);
- free(buff);
-
- uint_t dist_len = 4 * pow(Ls_c[i], 2);
- uint32_t *dist = malloc(dist_len * sizeof(uint32_t));
- FILE *dist_file = fopen(fn, "rb");
- fread(dist, sizeof(uint32_t), dist_len, dist_file);
- fclose(dist_file);
- {
- double mom2[2];
- mom(dist_len, 2, dist, mom2);
- vals_c2[i] = mom2[0];
- weights_c2[i] = mom2[1];
-
- double mom3[2];
- mom(dist_len, 3, dist, mom3);
- vals_c3[i] = mom3[0];
- weights_c3[i] = mom3[1];
- }
- free(dist);
- }
-
- const gsl_multifit_nlinear_type *T = gsl_multifit_nlinear_trust;
- gsl_multifit_nlinear_workspace *w;
- gsl_multifit_nlinear_fdf fdf;
- gsl_multifit_nlinear_parameters fdf_params = gsl_multifit_nlinear_default_parameters();
- uint_t n = 2 * nc;
- uint_t p = 27;
-
- double x_init[27] = { 1.56, .3, 2, -6, -10, -10, -10, 10, 15, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
- double weights[n];
- for (uint_t i = 0; i < nc; i++) {
- weights[i] = 1/gsl_pow_2(weights_c2[i] / vals_c2[i]);
- weights[nc + i] = 1/gsl_pow_2(weights_c3[i] / vals_c3[i]);
- }
-
- gsl_vector_view x = gsl_vector_view_array(x_init, p);
- gsl_vector_view wts = gsl_vector_view_array(weights, n);
-
- gsl_vector *f;
-
- fdf.f = f_moms;
- fdf.df = NULL;
- fdf.fvv = NULL;
- fdf.n = n;
- fdf.p = p;
- fdf.params = d;
-
- fdf_params.trs = gsl_multifit_nlinear_trs_lm;
-
- w = gsl_multifit_nlinear_alloc(T, &fdf_params, n, p);
- gsl_multifit_nlinear_winit(&x.vector, &wts.vector, &fdf, w);
- f = gsl_multifit_nlinear_residual(w);
- gsl_blas_ddot(f, f, &chisq0);
-
- const double xtol = 0.0;
- const double gtol = 0.0;
- const double ftol = 0.0;
-
- int info;
- int status = gsl_multifit_nlinear_driver(100, xtol, gtol, ftol, callback, NULL, &info, w);
- gsl_blas_ddot(f, f, &chisq);
-
- fprintf(stderr, "summary from method '%s/%s'\n",
- gsl_multifit_nlinear_name(w),
- gsl_multifit_nlinear_trs_name(w));
-
- fprintf(stderr, "number of iterations: %zu\n",
- gsl_multifit_nlinear_niter(w));
- fprintf(stderr, "function evaluations: %zu\n", fdf.nevalf);
- fprintf(stderr, "reason for stopping: %s\n",
- (info == 1) ? "small step size" : "small gradient");
- fprintf(stderr, "initial |f(x)| = %g\n", sqrt(chisq0));
- fprintf(stderr, "final |f(x)| = %g\n", sqrt(chisq));
-
- free(Ls_c);
- free(betas_c);
- free(vals_c2);
- free(weights_c2);
- free(vals_c3);
- free(weights_c3);
-
- return 0;
-}
-
diff --git a/src/fortune/defs.h b/src/fortune/defs.h
deleted file mode 100644
index 9f1b1db..0000000
--- a/src/fortune/defs.h
+++ /dev/null
@@ -1,134 +0,0 @@
-#ifndef NULL
-#define NULL 0
-#endif
-#define DELETED -2
-#include <stdlib.h>
-#include <limits.h>
-#include <stdbool.h>
-
-extern int triangulate, sorted, plot, debug;
-
-struct Freenode {
- struct Freenode *nextfree;
-};
-struct Freelist {
- struct Freenode *head;
- int nodesize;
-};
-char *getfree();
-char *myalloc();
-
-extern double xmin, xmax, ymin, ymax, deltax, deltay;
-
-struct Point {
- double x, y;
-};
-
-/* structure used both for sites and for vertices */
-struct Site {
- struct Point coord;
- int sitenbr;
- int refcnt;
-};
-
-extern struct Site *sites;
-extern int nsites;
-extern int siteidx;
-extern int sqrt_nsites;
-extern int nvertices;
-extern struct Freelist sfl;
-extern struct Site *bottomsite;
-
-void **alloclist;
-int allocnum;
-
-struct Edge {
- double a, b, c;
- struct Site *ep[2];
- struct Site *reg[2];
- int edgenbr;
-};
-#define le 0
-#define re 1
-extern int nedges;
-extern struct Freelist efl;
-
-int has_endpoint(), right_of();
-struct Site *intersect();
-double dist();
-struct Point PQ_min();
-struct Halfedge *PQextractmin();
-struct Edge *bisect();
-
-struct Halfedge {
- struct Halfedge *ELleft, *ELright;
- struct Edge *ELedge;
- int ELrefcnt;
- char ELpm;
- struct Site *vertex;
- double ystar;
- struct Halfedge *PQnext;
-};
-
-extern struct Freelist hfl;
-extern struct Halfedge *ELleftend, *ELrightend;
-extern int ELhashsize;
-extern struct Halfedge **ELhash;
-struct Halfedge *HEcreate(), *ELleft(), *ELright(), *ELleftbnd();
-struct Site *leftreg(), *rightreg();
-
-extern int PQhashsize;
-extern struct Halfedge *PQhash;
-extern struct Halfedge *PQfind();
-extern int PQcount;
-extern int PQmin;
-int PQempty();
-
-extern double *site_list;
-extern double *vert_list;
-extern double *line_list;
-extern unsigned int *edge_list;
-extern unsigned int *dual_list;
-
-extern int vert_count;
-extern int edge_count;
-extern int line_count;
-extern int dual_count;
-extern int site_count;
-
-void makefree(struct Freenode *curr, struct Freelist *fl);
-void voronoi(int triangulate, struct Site *(*nextsite)());
-char *getfree(struct Freelist *fl);
-char *myalloc(unsigned n);
-void freeinit(struct Freelist *fl, int size);
-void ELinitialize();
-struct Halfedge *HEcreate(struct Edge *e, int pm);
-void ELinsert(struct Halfedge *lb, struct Halfedge *new);
-struct Halfedge *ELgethash(int b);
-struct Halfedge *ELleftbnd(struct Point *p);
-void ELdelete(struct Halfedge *he);
-struct Halfedge *ELright(struct Halfedge *he);
-struct Halfedge *ELleft(struct Halfedge *he);
-struct Site *leftreg(struct Halfedge *he);
-struct Site *rightreg(struct Halfedge *he);
-void geominit();
-struct Edge *bisect(struct Site *s1, struct Site *s2);
-struct Site *intersect(struct Halfedge *el1, struct Halfedge *el2);
-int right_of(struct Halfedge *el, struct Point *p);
-void endpoint(struct Edge *e, int lr, struct Site *s);
-double dist(struct Site *s, struct Site *t);
-void makevertex(struct Site *v);
-void deref(struct Site *v);
-void ref(struct Site *v);
-void PQinsert(struct Halfedge *he, struct Site *v, double offset);
-void PQdelete(struct Halfedge *he);
-int PQbucket(struct Halfedge *he);
-int PQempty();
-struct Point PQ_min();
-struct Halfedge *PQextractmin();
-void PQinitialize();
-void out_bisector(struct Edge *e);
-void out_ep(struct Edge *e);
-void out_vertex(struct Site *v);
-void out_site(struct Site *s);
-void out_triple(struct Site *s1, struct Site *s2, struct Site *s3);
diff --git a/src/fortune/edgelist.c b/src/fortune/edgelist.c
deleted file mode 100644
index dfd85a7..0000000
--- a/src/fortune/edgelist.c
+++ /dev/null
@@ -1,136 +0,0 @@
-#
-#include "defs.h"
-int nedges;
-struct Freelist efl;
-struct Freelist hfl;
-struct Halfedge *ELleftend, *ELrightend;
-int ELhashsize;
-struct Halfedge **ELhash;
-
-int ntry, totalsearch;
-
-void ELinitialize() {
- int i;
- freeinit(&hfl, sizeof **ELhash);
- ELhashsize = 2 * sqrt_nsites;
- ELhash = (struct Halfedge **)myalloc(sizeof *ELhash * ELhashsize);
- for (i = 0; i < ELhashsize; i += 1)
- ELhash[i] = (struct Halfedge *)NULL;
- ELleftend = HEcreate((struct Edge *)NULL, 0);
- ELrightend = HEcreate((struct Edge *)NULL, 0);
- ELleftend->ELleft = (struct Halfedge *)NULL;
- ELleftend->ELright = ELrightend;
- ELrightend->ELleft = ELleftend;
- ELrightend->ELright = (struct Halfedge *)NULL;
- ELhash[0] = ELleftend;
- ELhash[ELhashsize - 1] = ELrightend;
-}
-
-struct Halfedge *HEcreate(e, pm) struct Edge *e;
-int pm;
-{
- struct Halfedge *answer;
- answer = (struct Halfedge *)getfree(&hfl);
- answer->ELedge = e;
- answer->ELpm = pm;
- answer->PQnext = (struct Halfedge *)NULL;
- answer->vertex = (struct Site *)NULL;
- answer->ELrefcnt = 0;
- return (answer);
-}
-
-void ELinsert(struct Halfedge *lb, struct Halfedge *new) {
- new->ELleft = lb;
- new->ELright = lb->ELright;
- (lb->ELright)->ELleft = new;
- lb->ELright = new;
-}
-
-/* Get entry from hash table, pruning any deleted nodes */
-struct Halfedge *ELgethash(b) int b;
-{
- struct Halfedge *he;
-
- if (b < 0 || b >= ELhashsize)
- return ((struct Halfedge *)NULL);
- he = ELhash[b];
- if (he == (struct Halfedge *)NULL || he->ELedge != (struct Edge *)DELETED)
- return (he);
-
- /* Hash table points to deleted half edge. Patch as necessary. */
- ELhash[b] = (struct Halfedge *)NULL;
- if ((he->ELrefcnt -= 1) == 0)
- makefree((struct Freenode *)he, &hfl);
- return ((struct Halfedge *)NULL);
-}
-
-struct Halfedge *ELleftbnd(p) struct Point *p;
-{
- int i, bucket;
- struct Halfedge *he;
-
- /* Use hash table to get close to desired halfedge */
- bucket = (p->x - xmin) / deltax * ELhashsize;
- if (bucket < 0)
- bucket = 0;
- if (bucket >= ELhashsize)
- bucket = ELhashsize - 1;
- he = ELgethash(bucket);
- if (he == (struct Halfedge *)NULL) {
- for (i = 1; 1; i += 1) {
- if ((he = ELgethash(bucket - i)) != (struct Halfedge *)NULL)
- break;
- if ((he = ELgethash(bucket + i)) != (struct Halfedge *)NULL)
- break;
- };
- totalsearch += i;
- };
- ntry += 1;
- /* Now search linear list of halfedges for the corect one */
- if (he == ELleftend || (he != ELrightend && right_of(he, p))) {
- do {
- he = he->ELright;
- } while (he != ELrightend && right_of(he, p));
- he = he->ELleft;
- } else
- do {
- he = he->ELleft;
- } while (he != ELleftend && !right_of(he, p));
-
- /* Update hash table and reference counts */
- if (bucket > 0 && bucket < ELhashsize - 1) {
- if (ELhash[bucket] != (struct Halfedge *)NULL)
- ELhash[bucket]->ELrefcnt -= 1;
- ELhash[bucket] = he;
- ELhash[bucket]->ELrefcnt += 1;
- };
- return (he);
-}
-
-/* This delete routine can't reclaim node, since pointers from hash
- table may be present. */
-void ELdelete(struct Halfedge *he) {
- (he->ELleft)->ELright = he->ELright;
- (he->ELright)->ELleft = he->ELleft;
- he->ELedge = (struct Edge *)DELETED;
-}
-
-struct Halfedge *ELright(he) struct Halfedge *he;
-{ return (he->ELright); }
-
-struct Halfedge *ELleft(he) struct Halfedge *he;
-{ return (he->ELleft); }
-
-struct Site *leftreg(he) struct Halfedge *he;
-{
- if (he->ELedge == (struct Edge *)NULL)
- return (bottomsite);
- return (he->ELpm == le ? he->ELedge->reg[le] : he->ELedge->reg[re]);
-}
-
-struct Site *rightreg(he) struct Halfedge *he;
-{
- if (he->ELedge == (struct Edge *)NULL)
- return (bottomsite);
- return (he->ELpm == le ? he->ELedge->reg[re] : he->ELedge->reg[le]);
-}
diff --git a/src/fortune/geometry.c b/src/fortune/geometry.c
deleted file mode 100644
index 131cacf..0000000
--- a/src/fortune/geometry.c
+++ /dev/null
@@ -1,184 +0,0 @@
-#
-#include "defs.h"
-#include <math.h>
-
-void geominit() {
- struct Edge e;
- double sn;
-
- freeinit(&efl, sizeof e);
- nvertices = 0;
- nedges = 0;
- sn = nsites + 4;
- sqrt_nsites = sqrt(sn);
- deltay = ymax - ymin;
- deltax = xmax - xmin;
-}
-
-struct Edge *bisect(s1, s2) struct Site *s1, *s2;
-{
- double dx, dy, adx, ady;
- struct Edge *newedge;
-
- newedge = (struct Edge *)getfree(&efl);
-
- newedge->reg[0] = s1;
- newedge->reg[1] = s2;
- ref(s1);
- ref(s2);
- newedge->ep[0] = (struct Site *)NULL;
- newedge->ep[1] = (struct Site *)NULL;
-
- dx = s2->coord.x - s1->coord.x;
- dy = s2->coord.y - s1->coord.y;
- adx = dx > 0 ? dx : -dx;
- ady = dy > 0 ? dy : -dy;
- newedge->c = s1->coord.x * dx + s1->coord.y * dy + (dx * dx + dy * dy) * 0.5;
- if (adx > ady) {
- newedge->a = 1.0;
- newedge->b = dy / dx;
- newedge->c /= dx;
- } else {
- newedge->b = 1.0;
- newedge->a = dx / dy;
- newedge->c /= dy;
- };
-
- newedge->edgenbr = nedges;
- out_bisector(newedge);
- nedges += 1;
- return (newedge);
-}
-
-struct Site *intersect(el1, el2) struct Halfedge *el1, *el2;
-{
- struct Edge *e1, *e2, *e;
- struct Halfedge *el;
- double d, xint, yint;
- int right_of_site;
- struct Site *v;
-
- e1 = el1->ELedge;
- e2 = el2->ELedge;
- if (e1 == (struct Edge *)NULL || e2 == (struct Edge *)NULL)
- return ((struct Site *)NULL);
- if (e1->reg[1] == e2->reg[1])
- return ((struct Site *)NULL);
-
- d = e1->a * e2->b - e1->b * e2->a;
- /* printf("intersect: d=%g\n", d); */
- if (-1.0e-20 < d && d < 1.0e-20) {
- return ((struct Site *)NULL);
- };
-
- xint = (e1->c * e2->b - e2->c * e1->b) / d;
- yint = (e2->c * e1->a - e1->c * e2->a) / d;
-
- if ((e1->reg[1]->coord.y < e2->reg[1]->coord.y) ||
- (e1->reg[1]->coord.y == e2->reg[1]->coord.y &&
- e1->reg[1]->coord.x < e2->reg[1]->coord.x)) {
- el = el1;
- e = e1;
- } else {
- el = el2;
- e = e2;
- };
- right_of_site = xint >= e->reg[1]->coord.x;
- if ((right_of_site && el->ELpm == le) || (!right_of_site && el->ELpm == re))
- return ((struct Site *)NULL);
-
- v = (struct Site *)getfree(&sfl);
- v->refcnt = 0;
- v->coord.x = xint;
- v->coord.y = yint;
- return (v);
-}
-
-/* returns 1 if p is to right of halfedge e */
-int right_of(el, p) struct Halfedge *el;
-struct Point *p;
-{
- struct Edge *e;
- struct Site *topsite;
- int right_of_site, above, fast;
- double dxp, dyp, dxs, t1, t2, t3, yl;
-
- e = el->ELedge;
- topsite = e->reg[1];
- right_of_site = p->x > topsite->coord.x;
- if (right_of_site && el->ELpm == le)
- return (1);
- if (!right_of_site && el->ELpm == re)
- return (0);
-
- if (e->a == 1.0) {
- dyp = p->y - topsite->coord.y;
- dxp = p->x - topsite->coord.x;
- fast = 0;
- if ((!right_of_site && e->b < 0.0) | (right_of_site && e->b >= 0.0)) {
- above = dyp >= e->b * dxp;
- fast = above;
- } else {
- above = p->x + p->y * e->b > e->c;
- if (e->b < 0.0)
- above = !above;
- if (!above)
- fast = 1;
- };
- if (!fast) {
- dxs = topsite->coord.x - (e->reg[0])->coord.x;
- above = e->b * (dxp * dxp - dyp * dyp) <
- dxs * dyp * (1.0 + 2.0 * dxp / dxs + e->b * e->b);
- if (e->b < 0.0)
- above = !above;
- };
- } else /*e->b==1.0 */
- {
- yl = e->c - e->a * p->x;
- t1 = p->y - yl;
- t2 = p->x - topsite->coord.x;
- t3 = yl - topsite->coord.y;
- above = t1 * t1 > t2 * t2 + t3 * t3;
- };
- return (el->ELpm == le ? above : !above);
-}
-
-void endpoint(e, lr, s) struct Edge *e;
-int lr;
-struct Site *s;
-{
- e->ep[lr] = s;
- ref(s);
- if (e->ep[re - lr] == (struct Site *)NULL) {
- } else {
- out_ep(e);
- deref(e->reg[le]);
- deref(e->reg[re]);
- makefree((struct Freenode *)e, &efl);
- }
-}
-
-double dist(s, t) struct Site *s, *t;
-{
- double dx, dy;
- dx = s->coord.x - t->coord.x;
- dy = s->coord.y - t->coord.y;
- return (sqrt(dx * dx + dy * dy));
-}
-
-void makevertex(v) struct Site *v;
-{
- v->sitenbr = nvertices;
- nvertices += 1;
- out_vertex(v);
-}
-
-void deref(v) struct Site *v;
-{
- v->refcnt -= 1;
- if (v->refcnt == 0)
- makefree((struct Freenode *)v, &sfl);
-}
-
-void ref(v) struct Site *v;
-{ v->refcnt += 1; }
diff --git a/src/fortune/heap.c b/src/fortune/heap.c
deleted file mode 100644
index 2ebbd58..0000000
--- a/src/fortune/heap.c
+++ /dev/null
@@ -1,94 +0,0 @@
-#
-#include "defs.h"
-int PQhashsize;
-struct Halfedge *PQhash;
-struct Halfedge *PQfind();
-int PQcount;
-int PQmin;
-int PQempty();
-
-void PQinsert(he, v, offset) struct Halfedge *he;
-struct Site *v;
-double offset;
-{
- struct Halfedge *last, *next;
- he->vertex = v;
- ref(v);
- he->ystar = v->coord.y + offset;
- last = &PQhash[PQbucket(he)];
- while ((next = last->PQnext) != (struct Halfedge *)NULL &&
- (he->ystar > next->ystar ||
- (he->ystar == next->ystar && v->coord.x > next->vertex->coord.x))) {
- last = next;
- };
- he->PQnext = last->PQnext;
- last->PQnext = he;
- PQcount += 1;
-}
-
-void PQdelete(he) struct Halfedge *he;
-{
- struct Halfedge *last;
-
- if (he->vertex != (struct Site *)NULL) {
- last = &PQhash[PQbucket(he)];
- while (last->PQnext != he)
- last = last->PQnext;
- last->PQnext = he->PQnext;
- PQcount -= 1;
- deref(he->vertex);
- he->vertex = (struct Site *)NULL;
- };
-}
-
-int PQbucket(he) struct Halfedge *he;
-{
- int bucket;
-
- if (he->ystar < ymin)
- bucket = 0;
- else if (he->ystar >= ymax)
- bucket = PQhashsize - 1;
- else
- bucket = (he->ystar - ymin) / deltay * PQhashsize;
- if (bucket < 0)
- bucket = 0;
- if (bucket >= PQhashsize)
- bucket = PQhashsize - 1;
- if (bucket < PQmin)
- PQmin = bucket;
- return (bucket);
-}
-
-int PQempty() { return (PQcount == 0); }
-
-struct Point PQ_min() {
- struct Point answer;
-
- while (PQhash[PQmin].PQnext == (struct Halfedge *)NULL) {
- PQmin += 1;
- };
- answer.x = PQhash[PQmin].PQnext->vertex->coord.x;
- answer.y = PQhash[PQmin].PQnext->ystar;
- return (answer);
-}
-
-struct Halfedge *PQextractmin() {
- struct Halfedge *curr;
- curr = PQhash[PQmin].PQnext;
- PQhash[PQmin].PQnext = curr->PQnext;
- PQcount -= 1;
- return (curr);
-}
-
-void PQinitialize() {
- int i;
-// struct Point *s;
-
- PQcount = 0;
- PQmin = 0;
- PQhashsize = 4 * sqrt_nsites;
- PQhash = (struct Halfedge *)myalloc(PQhashsize * sizeof *PQhash);
- for (i = 0; i < PQhashsize; i += 1)
- PQhash[i].PQnext = (struct Halfedge *)NULL;
-}
diff --git a/src/fortune/main.c b/src/fortune/main.c
deleted file mode 100644
index b6b4476..0000000
--- a/src/fortune/main.c
+++ /dev/null
@@ -1,287 +0,0 @@
-#
-#include <stdio.h>
-#include <stdint.h>
-#include "defs.h"
-#include <math.h>
-struct Site *nextone();
-
-int triangulate, sorted, plot, debug;
-double xmin, xmax, ymin, ymax, deltax, deltay;
-int vert_count, edge_count, line_count, dual_count, site_count;
-double *site_list, *vert_list, *line_list;
-unsigned int *edge_list, *dual_list;
-
-/* sort sites on y, then x, coord */
-int scomp(s1, s2) struct Point *s1, *s2;
-{
- if (s1->y < s2->y)
- return (-1);
- if (s1->y > s2->y)
- return (1);
- if (s1->x < s2->x)
- return (-1);
- if (s1->x > s2->x)
- return (1);
- return (0);
-}
-
-/* read all sites, sort, and compute xmin, xmax, ymin, ymax */
-void readsites(double *lattice) {
-
- sites = (struct Site *)myalloc(nsites * sizeof *sites);
- for (unsigned int j = 0; j < nsites; j++) {
- sites[j].coord.x = lattice[2 * j];
- sites[j].coord.y = lattice[2 * j + 1];
- sites[j].sitenbr = j;
- sites[j].refcnt = 0;
- }
- qsort(sites, nsites, sizeof(struct Site), scomp);
- xmin = sites[0].coord.x;
- xmax = sites[0].coord.x;
- for (unsigned int i = 1; i < nsites; i += 1) {
- if (sites[i].coord.x < xmin)
- xmin = sites[i].coord.x;
- if (sites[i].coord.x > xmax)
- xmax = sites[i].coord.x;
- };
- ymin = sites[0].coord.y;
- ymax = sites[nsites - 1].coord.y;
-}
-
-unsigned int delete_duplicates(unsigned int ne, unsigned int *etv) {
- unsigned int ndup = 0;
- bool *duplicates = (bool *)calloc(ne, sizeof(bool));
- for (unsigned int i = 0; i < ne; i++) {
- unsigned int v1 = etv[2 * i];
- unsigned int v2 = etv[2 * i + 1];
- for (unsigned int j = 0; j < i; j++) {
- unsigned int w1 = etv[2 * j];
- unsigned int w2 = etv[2 * j + 1];
- if (w1 == v1 && w2 == v2) {
- duplicates[j] = true;
- ndup++;
- break;
- }
- }
- }
- unsigned int newsize = (int)ne - (int)ndup;
- unsigned int count = 0;
- for (unsigned int i = 0; i < ne; i++) {
- if (!duplicates[i]) {
- etv[2 * count] = etv[2 * i];
- etv[2 * count + 1] = etv[2 * i + 1];
- count++;
- }
- }
- free(duplicates);
- return newsize;
-}
-
-intptr_t *run_voronoi(unsigned int num, double *lattice, bool periodic, double xmin, double xmax, double ymin, double ymax)
-{
- struct Site *(*next)();
-
- sorted = 0;
- triangulate = 0;
- plot = 0;
- debug = 0;
-
- alloclist = (void **)malloc(9 * num * sizeof(void *));
- allocnum = 0;
-
- freeinit(&sfl, sizeof *sites);
-
- unsigned int eff_num = num;
- double *eff_lattice = lattice;
-
- if (periodic) {
- eff_num = 9 * num;
- eff_lattice = (double *)malloc(2 * eff_num * sizeof(double));
-
- for (unsigned int i = 0; i < num; i++) {
- // original sites - our baby boys
- eff_lattice[2*i] = lattice[2*i];
- eff_lattice[2*i+1] = lattice[2*i+1];
-
- // sites to the right
- eff_lattice[2*(num+i)+1] = lattice[2*i+1] - xmin + xmax;
- eff_lattice[2*(num+i)] = lattice[2*i];
-
- // sites to the left
- eff_lattice[2*(2*num+i)+1] = lattice[2*i+1] - xmax + xmin;
- eff_lattice[2*(2*num+i)] = lattice[2*i];
-
- // sites to the top
- eff_lattice[2*(3*num+i)+1] = lattice[2*i+1];
- eff_lattice[2*(3*num+i)] = lattice[2*i] - ymin + ymax;
-
- // sites to the bottom
- eff_lattice[2*(4*num+i)+1] = lattice[2*i+1];
- eff_lattice[2*(4*num+i)] = lattice[2*i] - ymax + ymin;
-
- // sites to the upper right
- eff_lattice[2*(5*num+i)+1] = lattice[2*i+1] - xmin + xmax;
- eff_lattice[2*(5*num+i)] = lattice[2*i] - ymin + ymax;
-
- // sites to the upper left
- eff_lattice[2*(6*num+i)+1] = lattice[2*i+1] - xmax + xmin;
- eff_lattice[2*(6*num+i)] = lattice[2*i] - ymin + ymax;
-
- // sites to the lower left
- eff_lattice[2*(7*num+i)+1] = lattice[2*i+1] - xmax + xmin;
- eff_lattice[2*(7*num+i)] = lattice[2*i] + ymin - ymax;
-
- // sites to the lower right
- eff_lattice[2*(8*num+i)+1] = lattice[2*i+1] + xmax - xmin;
- eff_lattice[2*(8*num+i)] = lattice[2*i] + ymin - ymax;
- }
- }
-
- nsites = eff_num;
- readsites(eff_lattice);
- if (periodic) free(eff_lattice);
- next = nextone;
-
- siteidx = 0;
- geominit();
-
- site_list = malloc(2 * nsites * sizeof(double));
- vert_list = NULL;
- line_list = NULL;
- edge_list = NULL;
- dual_list = NULL;
- site_count = 0;
- vert_count = 0;
- edge_count = 0;
- line_count = 0;
- dual_count = 0;
-
- voronoi(triangulate, next);
-
- unsigned int real_vert_count = vert_count;
- unsigned int real_edge_count = edge_count;
- double *real_vert_list = vert_list;
- unsigned int *real_edge_list = edge_list;
- unsigned int *real_dual_list = dual_list;
-
- if (periodic) {
- real_vert_count = 0;
- real_edge_count = 0;
- real_vert_list = (double *)malloc(2 * vert_count * sizeof(double));
- real_edge_list = (unsigned int *)malloc(2 * edge_count * sizeof(unsigned int));
- real_dual_list = (unsigned int *)malloc(3 * dual_count * sizeof(unsigned int));
- unsigned int *triangle_analyzed = (unsigned int *)malloc(4 * dual_count * sizeof(unsigned int));
- unsigned int q = 0;
- int *new_vert_ind = (int *)malloc(vert_count * sizeof(int));
-
- for (unsigned int i = 0; i < vert_count; i++) {
- unsigned int t1, t2, t3;
- t1 = dual_list[3*i]; t2 = dual_list[3*i+1]; t3 = dual_list[3*i+2];
- if (t1 >= num && t2 >= num && t3 >= num) {
- new_vert_ind[i] = -1;
- }
- else if (t1 < num && t2 < num && t3 < num) {
- real_vert_list[2*real_vert_count] = vert_list[2*i];
- real_vert_list[2*real_vert_count+1] = vert_list[2*i+1];
- real_dual_list[3*real_vert_count] = t1;
- real_dual_list[3*real_vert_count+1] = t2;
- real_dual_list[3*real_vert_count+2] = t3;
- new_vert_ind[i] = real_vert_count;
- real_vert_count++;
- }
- else {
- unsigned int tt1, tt2, tt3;
- tt1 = t1 % num; tt2 = t2 % num; tt3 = t3 % num;
- unsigned int s1, s2, s3;
- s1 = tt1 < tt2 ? (tt1 < tt3 ? tt1 : tt3) : (tt2 < tt3 ? tt2 : tt3);
- s3 = tt1 > tt2 ? (tt1 > tt3 ? tt1 : tt3) : (tt2 > tt3 ? tt2 : tt3);
- s2 = tt1 < tt2 ? (tt1 > tt3 ? tt1 : (tt2 < tt3 ? tt2 : tt3)) : (tt1 < tt3 ? tt1 : (tt2 < tt3 ? tt3: tt2));
-
- bool matched = false;
- unsigned int ass_vert;
-
- for (unsigned int j = 0; j < q; j++) {
- if (s1 == triangle_analyzed[4*j] && s2 == triangle_analyzed[4*j+1] && s3 == triangle_analyzed[4*j+2]) {
- matched = true;
- ass_vert = triangle_analyzed[4*j+3];
- break;
- }
- }
-
- if (matched) {
- new_vert_ind[i] = ass_vert;
- } else {
- triangle_analyzed[4*q] = s1;
- triangle_analyzed[4*q+1] = s2;
- triangle_analyzed[4*q+2] = s3;
- triangle_analyzed[4*q+3] = real_vert_count;
- q++;
- real_vert_list[2*real_vert_count+1] = fmod(2*xmax+vert_list[2*i+1], xmax);
- real_vert_list[2*real_vert_count] = fmod(2*ymax+vert_list[2*i], ymax);
- real_dual_list[3*real_vert_count] = t1 % num;
- real_dual_list[3*real_vert_count+1] = t2 % num;
- real_dual_list[3*real_vert_count+2] = t3 % num;
- new_vert_ind[i] = real_vert_count;
- real_vert_count++;
- }
- }
- }
- for (unsigned int i = 0; i < edge_count; i++) {
- unsigned int v1, v2;
- v1 = edge_list[2 * i];
- v2 = edge_list[2 * i+1];
- if (v1 != UINT_MAX && v2 != UINT_MAX) {
- if (new_vert_ind[v1] >= 0 && new_vert_ind[v2] >=0) {
- unsigned int nv1 = new_vert_ind[v1];
- unsigned int nv2 = new_vert_ind[v2];
- real_edge_list[2*real_edge_count] = nv1 < nv2 ? nv1 : nv2;
- real_edge_list[2*real_edge_count+1] = nv1 < nv2 ? nv2 : nv1;
- real_edge_count++;
- }
- }
- }
- unsigned int new_edge_count = delete_duplicates(real_edge_count, real_edge_list);
- real_edge_count = new_edge_count;
- free(triangle_analyzed);
- free(new_vert_ind);
- free(dual_list);
- free(edge_list);
- free(vert_list);
- }
-
-
- intptr_t *output = (intptr_t *)malloc(6 * sizeof(intptr_t));
- output[0] = (intptr_t)malloc(4 * sizeof(unsigned int));
- ((unsigned int *)output[0])[0] = real_vert_count;
- ((unsigned int *)output[0])[1] = real_edge_count;
- ((unsigned int *)output[0])[2] = line_count;
- ((unsigned int *)output[0])[3] = real_vert_count;
- output[1] = (intptr_t)site_list;
- output[2] = (intptr_t)real_vert_list;
- output[3] = (intptr_t)real_edge_list;
- output[4] = (intptr_t)line_list;
- output[5] = (intptr_t)real_dual_list;
-
- free(ELhash);
- free(PQhash);
- free(sites);
- for (unsigned int i = 0; i < allocnum; i++) {
- free(alloclist[i]);
- }
- free(alloclist);
-
- return output;
-}
-
-/* return a single in-storage site */
-struct Site *nextone() {
- for (; siteidx < nsites; siteidx += 1) {
- if (siteidx == 0 || sites[siteidx].coord.x != sites[siteidx - 1].coord.x ||
- sites[siteidx].coord.y != sites[siteidx - 1].coord.y) {
- siteidx += 1;
- return (&sites[siteidx - 1]);
- };
- };
- return ((struct Site *)NULL);
-}
-
diff --git a/src/fortune/memory.c b/src/fortune/memory.c
deleted file mode 100644
index 3d62a92..0000000
--- a/src/fortune/memory.c
+++ /dev/null
@@ -1,44 +0,0 @@
-#
-#include "defs.h"
-#include <stdio.h>
-
-void freeinit(struct Freelist *fl, int size) {
- fl->head = (struct Freenode *)NULL;
- fl->nodesize = size;
-}
-
-char *getfree(fl) struct Freelist *fl;
-{
- int i;
- struct Freenode *t;
- if (fl->head == (struct Freenode *)NULL) {
- t = (struct Freenode *)myalloc(sqrt_nsites * fl->nodesize);
- alloclist[allocnum] = (void *)t;
- allocnum++;
- for (i = 0; i < sqrt_nsites; i += 1)
- makefree((struct Freenode *)((char *)t + i * fl->nodesize), fl);
- };
- t = fl->head;
- fl->head = (fl->head)->nextfree;
- return ((char *)t);
-}
-
-void makefree(curr, fl) struct Freenode *curr;struct Freelist *fl;
-{
- curr->nextfree = fl->head;
- fl->head = curr;
-}
-
-int total_alloc;
-char *myalloc(n) unsigned n;
-{
- char *t;
- if ((t = malloc(n)) == (char *)0) {
- fprintf(stderr,
- "Insufficient memory processing site %d (%d bytes in use)\n",
- siteidx, total_alloc);
- exit(1);
- };
- total_alloc += n;
- return (t);
-}
diff --git a/src/fortune/output.c b/src/fortune/output.c
deleted file mode 100644
index d496feb..0000000
--- a/src/fortune/output.c
+++ /dev/null
@@ -1,46 +0,0 @@
-#
-#include "defs.h"
-#include <stdio.h>
-double pxmin, pxmax, pymin, pymax, cradius;
-
-void out_bisector(e) struct Edge *e;
-{
- if (line_count % nsites == 0) line_list = realloc(line_list, 3 * (nsites + line_count) * sizeof(double));
- line_list[3*line_count] = e->a;
- line_list[3*line_count+1] = e->b;
- line_list[3*line_count+2] = e->c;
- line_count++;
-}
-
-void out_ep(e) struct Edge *e;
-{
- if (edge_count % nsites == 0) edge_list = realloc(edge_list, 2 * (nsites + edge_count) * sizeof(unsigned int));
- edge_list[2 * edge_count] = (e->ep[le] != (struct Site *)NULL ? e->ep[le]->sitenbr : UINT_MAX);
- edge_list[2 * edge_count + 1] = (e->ep[re] != (struct Site *)NULL ? e->ep[re]->sitenbr : UINT_MAX);
- edge_count++;
-}
-
-void out_vertex(v) struct Site *v;
-{
- if (vert_count % nsites == 0) vert_list = realloc(vert_list, 2 * (nsites + vert_count) * sizeof(double));
- vert_list[2 * vert_count] = v->coord.x;
- vert_list[2 * vert_count + 1] = v->coord.y;
- vert_count++;
-}
-
-void out_site(s) struct Site *s;
-{
- site_list[2 * site_count] = s->coord.x;
- site_list[2 * site_count + 1] = s->coord.y;
- site_count++;
-}
-
-void out_triple(s1, s2, s3) struct Site *s1, *s2, *s3;
-{
- if (dual_count % nsites == 0) dual_list = realloc(dual_list, 3 * (nsites + dual_count) * sizeof(unsigned int));
- dual_list[dual_count * 3] = s1->sitenbr;
- dual_list[dual_count * 3 + 1] = s2->sitenbr;
- dual_list[dual_count * 3 + 2] = s3->sitenbr;
- dual_count++;
-}
-
diff --git a/src/fortune/voronoi.c b/src/fortune/voronoi.c
deleted file mode 100644
index dc9945f..0000000
--- a/src/fortune/voronoi.c
+++ /dev/null
@@ -1,103 +0,0 @@
-#
-#include "defs.h"
-
-struct Site *sites;
-int nsites;
-int siteidx;
-int sqrt_nsites;
-int nvertices;
-struct Freelist sfl;
-struct Site *bottomsite;
-
-/* implicit parameters: nsites, sqrt_nsites, xmin, xmax, ymin, ymax,
- deltax, deltay (can all be estimates).
- Performance suffers if they are wrong; better to make nsites,
- deltax, and deltay too big than too small. (?) */
-
-void voronoi(triangulate, nextsite) int triangulate;
-struct Site *(*nextsite)();
-{
- struct Site *newsite, *bot, *top, *temp, *p;
- struct Site *v;
- struct Point newintstar;
- int pm;
- struct Halfedge *lbnd, *rbnd, *llbnd, *rrbnd, *bisector;
- struct Edge *e;
-
- PQinitialize();
- bottomsite = (*nextsite)();
- out_site(bottomsite);
- ELinitialize();
-
- newsite = (*nextsite)();
- while (1) {
- if (!PQempty())
- newintstar = PQ_min();
-
- if (newsite != (struct Site *)NULL &&
- (PQempty() || newsite->coord.y < newintstar.y ||
- (newsite->coord.y == newintstar.y &&
- newsite->coord.x < newintstar.x))) {/* new site is smallest */
- out_site(newsite);
- lbnd = ELleftbnd(&(newsite->coord));
- rbnd = ELright(lbnd);
- bot = rightreg(lbnd);
- e = bisect(bot, newsite);
- bisector = HEcreate(e, le);
- ELinsert(lbnd, bisector);
- if ((p = intersect(lbnd, bisector)) != (struct Site *)NULL) {
- PQdelete(lbnd);
- PQinsert(lbnd, p, dist(p, newsite));
- };
- lbnd = bisector;
- bisector = HEcreate(e, re);
- ELinsert(lbnd, bisector);
- if ((p = intersect(bisector, rbnd)) != (struct Site *)NULL) {
- PQinsert(bisector, p, dist(p, newsite));
- };
- newsite = (*nextsite)();
- } else if (!PQempty())
- /* intersection is smallest */
- {
- lbnd = PQextractmin();
- llbnd = ELleft(lbnd);
- rbnd = ELright(lbnd);
- rrbnd = ELright(rbnd);
- bot = leftreg(lbnd);
- top = rightreg(rbnd);
- out_triple(bot, top, rightreg(lbnd));
- v = lbnd->vertex;
- makevertex(v);
- endpoint(lbnd->ELedge, lbnd->ELpm, v);
- endpoint(rbnd->ELedge, rbnd->ELpm, v);
- ELdelete(lbnd);
- PQdelete(rbnd);
- ELdelete(rbnd);
- pm = le;
- if (bot->coord.y > top->coord.y) {
- temp = bot;
- bot = top;
- top = temp;
- pm = re;
- }
- e = bisect(bot, top);
- bisector = HEcreate(e, pm);
- ELinsert(llbnd, bisector);
- endpoint(e, re - pm, v);
- deref(v);
- if ((p = intersect(llbnd, bisector)) != (struct Site *)NULL) {
- PQdelete(llbnd);
- PQinsert(llbnd, p, dist(p, bot));
- };
- if ((p = intersect(bisector, rrbnd)) != (struct Site *)NULL) {
- PQinsert(bisector, p, dist(p, bot));
- };
- } else
- break;
- };
-
- for (lbnd = ELright(ELleftend); lbnd != ELrightend; lbnd = ELright(lbnd)) {
- e = lbnd->ELedge;
- out_ep(e);
- };
-}
diff --git a/src/fracture.c b/src/fracture.c
index 168e519..52db2d8 100644
--- a/src/fracture.c
+++ b/src/fracture.c
@@ -270,7 +270,7 @@ int main(int argc, char *argv[]) {
for (uint32_t i = 0; i < N; i++) {
printf("\033[F\033[JFRACTURE: %0*d / %d\n", (uint8_t)log10(N) + 1, i + 1, N);
- graph_t *g = graph_create(lattice, boundary, L, use_dual, &c);
+ graph_t *g = graph_create(lattice, boundary, L, use_dual);
net_t *net = net_create(g, inf, beta, crack_len, use_voltage_boundaries, &c);
net_t *tmp_net = net_copy(net, &c);
data_t *data = net_fracture(tmp_net, &c, cutoff);
@@ -406,7 +406,7 @@ int main(int argc, char *argv[]) {
}
if (save_cluster_dist) {
- uint_t *tmp_cluster_dist = get_cluster_dist(net, &c);
+ uint_t *tmp_cluster_dist = get_cluster_dist(net);
for (uint_t j = 0; j < g->dnv; j++) {
cluster_size_dist[j] += tmp_cluster_dist[j];
}
@@ -415,7 +415,7 @@ int main(int argc, char *argv[]) {
net_free(net, &c);
- graph_free(g, &c);
+ graph_free(g);
}
printf("\033[F\033[JFRACTURE: COMPLETE\n");
diff --git a/src/fracture.h b/src/fracture.h
index 0a3a687..b1114fb 100644
--- a/src/fracture.h
+++ b/src/fracture.h
@@ -23,6 +23,10 @@
#include <sys/types.h>
#include <unistd.h>
+#include <jst/graph.h>
+#include <jst/rand.h>
+
+
// these defs allow me to switch to long int cholmod in a sitch
#define int_t int
#define uint_t unsigned int
@@ -31,54 +35,6 @@
#define GSL_RAND_GEN gsl_rng_mt19937
-typedef enum lattice_t {
- VORONOI_LATTICE,
- SQUARE_LATTICE,
- VORONOI_HYPERUNIFORM_LATTICE
-} lattice_t;
-
-typedef enum bound_t {
- FREE_BOUND,
- CYLINDER_BOUND,
- TORUS_BOUND,
- EMBEDDED_BOUND
-} bound_t;
-
-typedef struct {
- uint_t ne;
- uint_t nv;
- uint_t dnv;
- uint_t *ev;
- uint_t *dev;
- double *vx;
- double *dvx;
-} frame_t;
-
-typedef struct {
- uint_t L;
- bound_t boundary;
- uint_t ne;
- uint_t nv;
- uint_t dnv;
- uint_t nb;
- uint_t *ev;
- uint_t *dev;
- double *vx;
- double *dvx;
- uint_t *vei;
- uint_t *ve;
- uint_t *dvei;
- uint_t *dve;
- uint_t *bi;
- uint_t *b;
- uint_t *nbi;
- uint_t *bni;
- bool *bq;
- uint_t *spanning_edges;
- uint_t num_spanning_edges;
- cholmod_sparse *voltcurmat;
-} graph_t;
-
typedef struct {
const graph_t *graph;
bool *fuses;
@@ -91,6 +47,7 @@ typedef struct {
uint_t dim;
uint_t nep;
uint_t *evp;
+ cholmod_sparse *voltcurmat;
} net_t;
typedef struct {
@@ -102,11 +59,6 @@ typedef struct {
intptr_t *run_voronoi(uint_t num_coords, double *coords, bool periodic, double xmin, double xmax, double ymin, double ymax);
-int update_components(const cholmod_sparse *laplacian, uint_t *marks,
- int old_num_components, int v1, int v2, int exclude);
-
-uint_t *find_components(const cholmod_sparse *laplacian, uint_t skip);
-
cholmod_sparse *gen_adjacency(const net_t *net, bool dual, bool use_gp, bool symmetric, cholmod_common *c);
cholmod_sparse *gen_laplacian(const net_t *net, cholmod_common *c);
@@ -121,6 +73,7 @@ double dual_vert_to_coord(uint_t width, bool periodic, uint_t vert,
bool index);
void factor_update(cholmod_factor *factor, uint_t v1, uint_t v2, cholmod_common *c);
+void factor_update2(cholmod_factor *factor, uint_t v1, cholmod_common *c);
void net_notch(net_t *net, double notch_len, cholmod_common *c);
data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff);
@@ -142,35 +95,19 @@ cholmod_sparse *gen_voltcurmat(uint_t num_edges, uint_t num_verts,
net_t *net_copy(const net_t *net, cholmod_common *c);
-graph_t *ini_square_network(uint_t width, bound_t boundary, bool side_bounds,
- cholmod_common *c);
-
-void graph_free(graph_t *network, cholmod_common *c);
void net_free(net_t *instance, cholmod_common *c);
net_t *net_create(const graph_t *g, double inf, double beta, double notch_len, bool vb, cholmod_common *c);
-graph_t *ini_voro_graph(uint_t L, bound_t boundary, bool use_dual,
- double *(*genfunc)(uint_t, bound_t, gsl_rng *, uint_t *),
- cholmod_common *c);
-
bool break_edge(net_t *instance, uint_t edge, cholmod_common *c);
-void finish_instance(net_t *instance, cholmod_common *c);
-
-net_t *coursegrain_square(net_t *instance, graph_t *network_p, cholmod_common *c);
-
-uint_t *get_clusters(net_t *instance, cholmod_common *c);
+components_t *get_clusters(net_t *instance);
-uint_t *get_cluster_dist(net_t *instance, cholmod_common *c);
+uint_t *get_cluster_dist(net_t *instance);
-double *genfunc_uniform(uint_t num, gsl_rng *r);
-double *genfunc_hyperuniform(uint_t num, gsl_rng *r);
void randfunc_flat(gsl_rng *r, double *x, double *y);
void randfunc_gaus(gsl_rng *r, double *x, double *y);
-uint_t **get_dists(const graph_t *network);
-uint_t *dijkstra(const graph_t *g, uint_t source);
double *get_corr(net_t *instance, uint_t **dists, cholmod_common *c);
double *bin_values(graph_t *network, uint_t width, double *values);
@@ -181,18 +118,6 @@ data_t *data_create(uint_t num_edges);
void data_free(data_t *data);
void data_update(data_t *data, uint_t last_broke, long double strength, double conductivity);
-graph_t *graph_create(lattice_t lattice, bound_t bound, uint_t L, bool dual, cholmod_common *c);
-
-uint_t find_cycles(uint_t num_edges, const bool *fuses, const uint_t *ev, const uint_t *vei, const uint_t *ve, int **cycles);
-
-bool set_connected(const cholmod_sparse *laplacian, uint_t *marks, int vertex, int label, int stop_at, int exclude);
-
-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);
-
-frame_t *frame_create(lattice_t lattice, uint_t L, bool dual);
-void frame_free(frame_t *frame);
bool is_in(uint_t len, uint_t *list, uint_t element);
diff --git a/src/frame_create.c b/src/frame_create.c
deleted file mode 100644
index 89ce69d..0000000
--- a/src/frame_create.c
+++ /dev/null
@@ -1,167 +0,0 @@
-
-#include "fracture.h"
-
-uint_t *get_voro_dual_edges(uint_t num_edges,
- uint_t num_verts, uint_t *edges,
- uint_t *triangles) {
- uint_t *dual_edges =
- (uint_t *)malloc(2 * num_edges * sizeof(uint_t));
- uint_t place = 0;
- for (uint_t i = 0; i < num_edges; i++) {
- uint_t v1, v2;
- v1 = edges[2 * i];
- v2 = edges[2 * i + 1];
- if (v1 < num_verts && v2 < num_verts) {
- bool found_match = false;
- for (uint_t j = 0; j < 3; j++) {
- for (uint_t k = 0; k < 3; k++) {
- uint_t t11, t12, t21, t22;
- t11 = triangles[3 * v1 + j];
- t12 = triangles[3 * v1 + ((j + 1) % 3)];
- t21 = triangles[3 * v2 + k];
- t22 = triangles[3 * v2 + ((k + 1) % 3)];
- if ((t11 == t21 && t12 == t22) || (t11 == t22 && t12 == t21)) {
- dual_edges[2 * place] = t11 < t12 ? t11 : t12;
- dual_edges[2 * place + 1] = t11 < t12 ? t12 : t11;
- place++;
- found_match = true;
- break;
- }
- }
- if (found_match)
- break;
- }
- }
- }
-
- return dual_edges;
-}
-
-frame_t *frame_create_voronoi(uint_t L, bool dual, bool hyperuniform) {
- double *dvx = NULL;
- uint_t dnv = 2 * pow(L / 2, 2);
-
- {
- gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
- gsl_rng_set(r, rand_seed());
-
- if (hyperuniform) {
- dvx = genfunc_hyperuniform(dnv, r);
- } else {
- dvx = genfunc_uniform(dnv, r);
- }
-
- gsl_rng_free(r);
- }
-
- // retrieve a periodic voronoi tesselation of the lattice
- intptr_t *vout = run_voronoi(dnv, dvx, true, 0, 1, 0, 1);
-
- uint_t nv = ((uint_t *)vout[0])[0];
- uint_t ne = ((uint_t *)vout[0])[1];
- double *vx = (double *)vout[2];
- uint_t *ev = (uint_t *)vout[3];
- uint_t *voro_tris = (uint_t *)vout[5];
-
- free((void *)vout[0]);
- free((void *)vout[1]);
- free((void *)vout[4]);
- free(vout);
-
- // get dual edges of the fully periodic graph
- uint_t *dev = get_voro_dual_edges(ne, nv, ev, voro_tris);
-
- frame_t *frame = (frame_t *)malloc(sizeof(frame_t));
- frame->ne = ne;
-
- // when use_dual is specificed, the edge and vertex sets are swapped with the
- // dual edge and dual vertex sets. once formally relabelled, everything
- // works the same way
- if (dual) {
- frame->nv = dnv;
- frame->dnv = nv;
- frame->ev = dev;
- frame->dev = ev;
- frame->vx = dvx;
- frame->dvx = vx;
- } else {
- frame->nv = nv;
- frame->dnv = dnv;
- frame->ev = ev;
- frame->dev = dev;
- frame->vx = vx;
- frame->dvx = dvx;
- }
-
- return frame;
-}
-
-frame_t *frame_create_square(uint_t L, bool dual) {
- uint_t nv = 2 * pow(L / 2, 2);
- uint_t ne = pow(L, 2);
-
- uint_t *ev = (uint_t *)malloc(2 * ne * sizeof(uint_t));
- uint_t *dev = (uint_t *)malloc(2 * ne * sizeof(uint_t));
- double *vx = (double *)malloc(2 * nv * sizeof(double));
- double *dvx = (double *)malloc(2 * nv * sizeof(double));
-
- for (uint_t i = 0; i < ne; i++) {
- uint_t x = i / L;
- uint_t y = i % L;
-
- ev[2 * i] = (L * x) / 2 + ((y + x % 2) / 2) % (L / 2);
- ev[2 * i + 1] = ((L * (x + 1)) / 2 + ((y + (x + 1) % 2) / 2) % (L / 2)) % nv;
-
- dev[2 * i] = (L * x) / 2 + ((y + (x + 1) % 2) / 2) % (L / 2);
- dev[2 * i + 1] = ((L * (x + 1)) / 2 + ((y + x % 2) / 2) % (L / 2)) % nv;
- }
-
- double dx = 1. / L;
-
- for (uint_t i = 0; i < nv; i++) {
- vx[2 * i] = dx * ((1 + i / (L / 2)) % 2 + 2 * (i % (L / 2)));
- vx[2 * i + 1] = dx * (i / (L / 2));
-
- dvx[2 * i] = dx * ((i / (L / 2)) % 2 + 2 * (i % (L / 2)));
- dvx[2 * i + 1] = dx * (i / (L / 2));
- }
-
- frame_t *frame = (frame_t *)malloc(sizeof(frame_t));
- frame->ne = ne;
- frame->nv = nv;
- frame->dnv = nv;
-
- if (dual) {
- frame->ev = dev;
- frame->dev = ev;
- frame->vx = dvx;
- frame->dvx = vx;
- } else {
- frame->ev = ev;
- frame->dev = dev;
- frame->vx = vx;
- frame->dvx = dvx;
- }
-
- return frame;
-}
-
-frame_t *frame_create(lattice_t lattice, uint_t L, bool dual) {
- switch (lattice) {
- case (SQUARE_LATTICE):
- return frame_create_square(L, dual);
- case (VORONOI_LATTICE):
- return frame_create_voronoi(L, dual, false);
- case (VORONOI_HYPERUNIFORM_LATTICE):
- return frame_create_voronoi(L, dual, true);
- }
-}
-
-void frame_free(frame_t *frame) {
- free(frame->ev);
- free(frame->dev);
- free(frame->vx);
- free(frame->dvx);
- free(frame);
-}
-
diff --git a/src/gen_laplacian.c b/src/gen_laplacian.c
deleted file mode 100644
index 1cc93a4..0000000
--- a/src/gen_laplacian.c
+++ /dev/null
@@ -1,137 +0,0 @@
-
-#include "fracture.h"
-
-cholmod_sparse *gen_adjacency(const net_t *net, bool dual, bool use_gp, bool symmetric, cholmod_common *c) {
- const graph_t *g = net->graph;
- uint_t nv;
- uint_t ne;
- uint_t nre;
- uint_t *ev;
-
- if (use_gp) {
- nv = net->dim;
- ne = net->nep;
- nre = (int_t)net->nep - (int_t)net->num_broken;
- ev = net->evp;
- } else if (dual) {
- nv = g->dnv;
- ne = g->ne;
- nre = net->num_broken;
- ev = g->dev;
- } else {
- nv = g->nv;
- ne = g->ne;
- nre = (int_t)g->ne - (int_t)net->num_broken;
- ev = g->ev;
- }
-
- uint_t nnz = nre;
-
- cholmod_triplet *t = CHOL_F(allocate_triplet)(nv, nv, nnz, 1, CHOLMOD_REAL, c);
-
- int_t *ri = (int_t *)t->i;
- int_t *ci = (int_t *)t->j;
- double *ai = (double *)t->x;
-
- t->nnz = nnz;
-
- uint_t a = 0;
-
- for (uint_t i = 0; i < ne; i++) {
- if ((net->fuses[i] && dual) || (!net->fuses[i] && !dual)) {
- uint_t v1 = ev[2 * i];
- uint_t v2 = ev[2 * i + 1];
-
- uint_t s1 = v1 < v2 ? v1 : v2;
- uint_t s2 = v1 < v2 ? v2 : v1;
-
- ri[a] = s2;
- ci[a] = s1;
- ai[a] = 1;
-
- a++;
- }
- }
-
- cholmod_sparse *s = CHOL_F(triplet_to_sparse)(t, nnz, c);
- CHOL_F(free_triplet)(&t, c);
-
- if (!symmetric) {
- cholmod_sparse *tmp_s = CHOL_F(copy)(s, 0, 1, c);
- CHOL_F(free_sparse)(&s, c);
- s = tmp_s;
- }
-
- return s;
-}
-
-cholmod_sparse *gen_laplacian(const net_t *net, cholmod_common *c) {
- const graph_t *g = net->graph;
- uint_t nv = net->dim;
- uint_t ne = net->nep;
- uint_t *ev = net->evp;
-
- uint_t nnz = nv;
-
- cholmod_triplet *temp_m = CHOL_F(allocate_triplet)(nv, nv, nnz, 1, CHOLMOD_REAL, c);
-
- int_t *rowind = (int_t *)temp_m->i;
- int_t *colind = (int_t *)temp_m->j;
- double *acoo = (double *)temp_m->x;
-
- temp_m->nnz = nnz;
-
- for (uint_t i = 0; i < nv; i++) {
- rowind[i] = i;
- colind[i] = i;
- acoo[i] = 0;
- }
-
- cholmod_sparse *adjacency = gen_adjacency(net, false, true, true, c);
-
- for (uint_t i = 0; i < ne; i++) {
- if (!net->fuses[i]){
- uint_t v1 = ev[2 * i];
- uint_t v2 = ev[2 * i + 1];
-
- acoo[v1]++;
- acoo[v2]++;
- }
- }
-
- if (net->voltage_bound && g->boundary != TORUS_BOUND) {
- for (uint_t i = 0; i < net->dim; i++) {
- uint_t v = g->nbi[i];
- for (uint_t j = 0; j < g->vei[v+1] - g->vei[v]; j++) {
- uint_t e = g->ve[g->vei[v] + j];
- uint_t v0 = g->ev[2 * e];
- uint_t v1 = g->ev[2 * e + 1];
-
- if (g->bq[v0] || g->bq[v1]) {
- acoo[i]++;
- }
- }
- }
- } else {
- acoo[0]++;
- }
-
- for (uint_t i = 0; i < nv; i++) {
- if (acoo[i] == 0) acoo[i]++;
- }
-
- //assert(CHOL_F(check_triplet)(temp_m, c));
-
- cholmod_sparse *t_out = CHOL_F(triplet_to_sparse)(temp_m, nnz, c);
- //assert(CHOL_F(check_sparse)(t_out, c));
-
- double alpha[2] = {1, 0};
- double beta[2] = {-1, 0};
- cholmod_sparse *laplacian = CHOL_F(add)(t_out, adjacency, alpha, beta, 1, 1, c);
-
- CHOL_F(free_sparse)(&t_out, c);
- CHOL_F(free_sparse)(&adjacency, c);
- CHOL_F(free_triplet)(&temp_m, c);
-
- return laplacian;
-}
diff --git a/src/gen_voltcurmat.c b/src/gen_voltcurmat.c
deleted file mode 100644
index e870140..0000000
--- a/src/gen_voltcurmat.c
+++ /dev/null
@@ -1,36 +0,0 @@
-
-#include "fracture.h"
-
-cholmod_sparse *gen_voltcurmat(unsigned int num_edges, unsigned int num_verts,
- unsigned int *edges_to_verts,
- cholmod_common *c) {
-
- cholmod_triplet *t_m = CHOL_F(allocate_triplet)(
- num_edges, num_verts, num_edges * 2, 0, CHOLMOD_REAL, c);
- assert(t_m != NULL);
-
- int_t *rowind = (int_t *)t_m->i;
- int_t *colind = (int_t *)t_m->j;
- double *acoo = (double *)t_m->x;
-
- for (unsigned int i = 0; i < num_edges; i++) {
- unsigned int v1 = edges_to_verts[2 * i];
- unsigned int v2 = edges_to_verts[2 * i + 1];
- rowind[2 * i] = i;
- rowind[2 * i + 1] = i;
- colind[2 * i] = v1;
- colind[2 * i + 1] = v2;
- acoo[2 * i] = 1;
- acoo[2 * i + 1] = -1;
- }
-
- t_m->nnz = num_edges * 2;
-
- assert(CHOL_F(check_triplet)(t_m, c));
-
- cholmod_sparse *m = CHOL_F(triplet_to_sparse)(t_m, num_edges * 2, c);
-
- CHOL_F(free_triplet)(&t_m, c);
-
- return m;
-}
diff --git a/src/geometry.c b/src/geometry.c
deleted file mode 100644
index ec788f1..0000000
--- a/src/geometry.c
+++ /dev/null
@@ -1,55 +0,0 @@
-
-#include "fracture.h"
-
-int edge_to_verts(unsigned int width, bool periodic, unsigned int edge,
- bool index) {
- assert(edge < pow(width, 2));
-
- int x = edge / width + 1;
- int y = edge % width + 1;
-
- if (periodic) {
- return (((index ^ (x % 2)) + 2 * ((y + (index ^ (!(x % 2)))) / 2) - 1) %
- width +
- (x - index) * width) /
- 2;
- } else {
- return ((index ^ (x % 2)) + 2 * ((y + (index ^ (!(x % 2)))) / 2) +
- (x - index) * (width + 1) - 1) /
- 2;
- }
-}
-
-int dual_edge_to_verts(unsigned int width, bool periodic, unsigned int edge,
- bool index) {
- assert(edge < pow(width, 2));
-
- int x = edge / width + 1;
- int y = edge % width + 1;
-
- if (periodic) {
- return (((index ^ (!(x % 2))) + 2 * ((y + (index ^ (x % 2))) / 2) - 1) %
- width +
- (x - index) * width) /
- 2;
- } else {
- return ((index ^ (!(x % 2))) + 2 * ((y + (index ^ (x % 2))) / 2) +
- (x - index) * (width + 1) - 1) /
- 2;
- }
-}
-
-double dual_vert_to_coord(unsigned int width, bool periodic, unsigned int vert,
- bool index) {
- if (periodic) {
- if (index)
- return (2 * vert) % width + (2 * vert / width) % 2;
- else
- return 2 * vert / width;
- } else {
- if (index)
- return (2 * vert) % (width + 1);
- else
- return (2 * vert) / (width + 1);
- }
-}
diff --git a/src/get_dual_clusters.c b/src/get_dual_clusters.c
deleted file mode 100644
index 6090fb8..0000000
--- a/src/get_dual_clusters.c
+++ /dev/null
@@ -1,36 +0,0 @@
-
-#include "fracture.h"
-
-unsigned int *get_clusters(net_t *instance, cholmod_common *c) {
- cholmod_sparse *s_dual = gen_adjacency(instance, true, false, false, c);
-
- unsigned int *dual_marks = find_components(s_dual, 0);
- CHOL_F(free_sparse)(&s_dual, c);
-
- return dual_marks;
-}
-
-unsigned int *get_cluster_dist(net_t *instance, cholmod_common *c) {
- unsigned int *clusters = get_clusters(instance, c);
- unsigned int *cluster_dist = (unsigned int *)calloc(
- instance->graph->dnv, sizeof(unsigned int));
-
- unsigned int cur_mark = 0;
- while (true) {
- cur_mark++;
- unsigned int num_in_cluster = 0;
- for (unsigned int i = 0; i < instance->graph->dnv; i++) {
- if (clusters[i] == cur_mark)
- num_in_cluster++;
- }
-
- if (num_in_cluster == 0)
- break;
-
- cluster_dist[num_in_cluster - 1]++;
- }
-
- free(clusters);
-
- return cluster_dist;
-}
diff --git a/src/graph_components.c b/src/graph_components.c
deleted file mode 100644
index 9680675..0000000
--- a/src/graph_components.c
+++ /dev/null
@@ -1,185 +0,0 @@
-
-#include "fracture.h"
-
-bool set_connected(const cholmod_sparse *laplacian, uint_t *marks,
- int vertex, int label, int stop_at, int exclude) {
- unsigned int num_verts = (int_t)laplacian->nrow;
-
- int_t *ia = (int_t *)laplacian->p;
- int_t *ja = (int_t *)laplacian->i;
- double *a = (double *)laplacian->x;
-
- bool stopped = false;
-
- int *queue = (int *)malloc((num_verts) * sizeof(int));
- assert(queue != NULL);
-
- int queue_pos = 0;
-
- int *component_list = (int *)malloc((num_verts) * sizeof(int));
- assert(component_list != NULL);
-
- int list_pos = 0;
- queue[0] = vertex;
- component_list[0] = vertex;
-
- unsigned int old_label = marks[vertex];
- marks[vertex] = label;
-
- while (queue_pos >= 0) {
- int next_vertex = queue[queue_pos];
- queue_pos--;
- for (int j = 0; j < ia[next_vertex + 1] - ia[next_vertex]; j++) {
- if (ja[ia[next_vertex] + j] == stop_at && a[ia[next_vertex] + j] != 0) {
- // if we run into the other vertex, the components haven't changed
- stopped = true; queue_pos = -1;
- break;
- }
- if (marks[ja[ia[next_vertex] + j]] != label &&
- a[ia[next_vertex] + j] != 0 &&
- ja[ia[next_vertex] + j] < num_verts - exclude) {
- marks[ja[ia[next_vertex] + j]] = label;
- queue_pos++;
- queue[queue_pos] = ja[ia[next_vertex] + j];
- list_pos++;
- component_list[list_pos] = ja[ia[next_vertex] + j];
- }
- }
- }
-
- if (stopped) {
- for (unsigned int i = 0; i <= list_pos; i++) {
- marks[component_list[i]] = old_label;
- }
- }
-
- free(queue);
- free(component_list);
-
- if (stopped)
- return false;
- else
- return true;
-}
-
-int update_components(const cholmod_sparse *laplacian, unsigned int *marks,
- int old_num_components, int v1, int v2, int exclude) {
- assert(laplacian != NULL);
- assert(marks != NULL);
-
- if (set_connected(laplacian, marks, v1, old_num_components + 1, v2,
- exclude)) {
- return old_num_components + 1;
- } else
- return old_num_components;
-}
-
-unsigned int *find_components(const cholmod_sparse *laplacian, unsigned int skip) {
- size_t num_verts = laplacian->nrow;
- unsigned int *marks = (unsigned int *)calloc(num_verts, sizeof(unsigned int));
-
- int num_components = 0;
-
- for (int i = 0; i < num_verts; i++) {
- if (marks[i] == 0) {
- num_components++;
- set_connected(laplacian, marks, i, num_components, -1, skip);
- }
- }
-
- return marks;
-}
-
-uint_t find_cycles(uint_t num_edges, const bool *fuses, const uint_t *ev, const unsigned
- int *vei, const uint_t *ve, int **cycles) {
- uint_t num_cycles = 0;
- uint_t *elist = (uint_t *)malloc(num_edges * sizeof(uint_t));
- uint_t *side = (uint_t *)calloc(num_edges, sizeof(uint_t));
- uint_t *num = (uint_t *)malloc(num_edges * sizeof(uint_t));
- for (uint_t i = 0; i < num_edges; i++) {
- if (fuses[i]) {
- bool in_cycle = false;
- for (uint_t j = 0; j < num_cycles; j++) {
- uint_t k = 0;
- int e;
- while ((e = cycles[j][k]) >= 0) {
- if (e == i) {
- in_cycle = true;
- break;
- }
- k++;
- }
- if (in_cycle) {
- break;
- }
- }
- if (!in_cycle) {
- uint_t pos = 0;
- uint_t cur_e = i;
- bool *included = (bool *)calloc(num_edges, sizeof(bool));
- included[cur_e] = true;
- elist[0] = cur_e;
- side[0] = 1;
- num[0] = 0;
- while (true) {
- uint_t cv = ev[2 * cur_e + side[pos]];
- uint_t new_e = cur_e;
- uint_t j;
- for (j = num[pos]; j < vei[cv+1] - vei[cv]; j++) {
- new_e = ve[vei[cv] + j];
- if (new_e != cur_e && fuses[new_e]) break;
- }
- if (new_e != cur_e && fuses[new_e]) {
- if (new_e == elist[0] && cv == ev[2 * elist[0] + !side[0]]) {
- pos++;
- cycles[2*num_cycles] = (int *)malloc((pos + 1) * sizeof(int));
- cycles[2*num_cycles+1] = (int *)malloc((pos + 1) * sizeof(int));
- for (uint_t i = 0; i < pos; i++) {
- cycles[2*num_cycles][i] = elist[i];
- cycles[2*num_cycles+1][i] = side[i];
- }
- cycles[2*num_cycles][pos] = -1;
- cycles[2*num_cycles+1][pos] = -1;
- num_cycles++;
- pos--;
- num[pos]++;
- included[pos] = false;
- cur_e = elist[pos];
- } else if (included[new_e]) {
- if (pos == 0) break;
- included[cur_e] = false;
- pos--;
- num[pos]++;
- cur_e = elist[pos];
- } else {
- num[pos] = j;
- pos++;
- elist[pos] = new_e;
- included[new_e] = true;
- uint_t nv1 = ev[2 * new_e];
- uint_t nv2 = ev[2 * new_e + 1];
- uint_t v = nv1 == cv ? nv2 : nv1;
- side[pos] = v == nv1 ? 0 : 1;
- num[pos] = 0;
- cur_e = new_e;
- }
- } else {
- if (pos == 0) break;
- included[cur_e] = false;
- pos--;
- num[pos]++;
- cur_e = elist[pos];
- }
- }
- free(included);
- }
- }
- }
-
- free(elist);
- free(side);
- free(num);
-
- return num_cycles;
-}
-
diff --git a/src/graph_create.c b/src/graph_create.c
deleted file mode 100644
index 09a3aed..0000000
--- a/src/graph_create.c
+++ /dev/null
@@ -1,299 +0,0 @@
-
-#include "fracture.h"
-
-uint_t *get_spanning_edges(uint_t num_edges, uint_t *edges_to_verts, double *vert_coords, double cut, uint_t *n) {
- uint_t *spanning_edges = (uint_t *)malloc(num_edges * sizeof(uint_t));
- (*n) = 0;
- for (uint_t i = 0; i < num_edges; i++) {
- uint_t v1, v2;
- v1 = edges_to_verts[2 * i];
- v2 = edges_to_verts[2 * i + 1];
- double v1y, v2y;
- v1y = vert_coords[2 * v1 + 1];
- v2y = vert_coords[2 * v2 + 1];
- if ((fabs(v1y - v2y) < 0.5) && ((v1y <= cut && v2y > cut) || (v1y >= cut && v2y < cut))) {
- spanning_edges[*n] = i;
- (*n)++;
- }
- }
- return spanning_edges;
-}
-
-double *get_edge_coords(uint_t num_edges, double *vert_coords,
- uint_t *edges_to_verts) {
- double *output = (double *)malloc(2 * num_edges * sizeof(double));
-
- for (uint_t i = 0; i < num_edges; i++) {
- uint_t v1, v2;
- double v1x, v1y, v2x, v2y, dx, dy;
- v1 = edges_to_verts[2 * i];
- v2 = edges_to_verts[2 * i + 1];
- output[2 * i] = 0;
- output[2 * i + 1] = 0;
- v1x = vert_coords[2 * v1];
- v1y = vert_coords[2 * v1 + 1];
- v2x = vert_coords[2 * v2];
- v2y = vert_coords[2 * v2 + 1];
- dx = v1x - v2x;
- dy = v1y - v2y;
- if (fabs(dx) > 0.5) {
- if (dx > 0) {
- v2x = v2x + 1;
- } else {
- v1x = v1x + 1;
- }
- }
- if (fabs(dy) > 0.5) {
- if (dy > 0) {
- v2y = v2y + 1;
- } else {
- v1y = v1y + 1;
- }
- }
- output[2 * i] = (v1x + v2x) / 2;
- output[2 * i + 1] = (v1y + v2y) / 2;
- }
-
- return output;
-}
-
-uint_t *get_verts_to_edges_ind(uint_t num_verts,
- uint_t num_edges,
- const uint_t *edges_to_verts) {
- uint_t *output =
- (uint_t *)calloc(num_verts + 1, sizeof(uint_t));
- assert(output != NULL);
-
- for (uint_t i = 0; i < 2 * num_edges; i++) {
- if (edges_to_verts[i] < num_verts) {
- output[edges_to_verts[i] + 1]++;
- }
- }
-
- for (uint_t i = 0; i < num_verts; i++) {
- output[i + 1] += output[i];
- }
-
- return output;
-}
-
-uint_t *get_verts_to_edges(uint_t num_verts, uint_t num_edges,
- const uint_t *edges_to_verts,
- const uint_t *verts_to_edges_ind) {
- uint_t *output = (uint_t *)calloc(verts_to_edges_ind[num_verts],
- sizeof(uint_t));
- uint_t *counts =
- (uint_t *)calloc(num_verts, sizeof(uint_t));
- for (int i = 0; i < 2 * num_edges; i++) {
- if (edges_to_verts[i] < num_verts) {
- output[verts_to_edges_ind[edges_to_verts[i]] +
- counts[edges_to_verts[i]]] = i / 2;
- counts[edges_to_verts[i]]++;
- }
- }
-
- free(counts);
-
- return output;
-}
-
-uint_t get_cut_edges(uint_t ne, const uint_t *ev, const double *vx, bool both, uint_t *ce) {
- uint_t nce = 0;
-
- for (uint_t i = 0; i < ne; i++) {
- uint_t v1 = ev[2 * i];
- uint_t v2 = ev[2 * i + 1];
-
- double v1y = vx[2 * v1 + 1];
- double v2y = vx[2 * v2 + 1];
-
- if (fabs(v1y - v2y) > 0.5) {
- ce[nce] = i;
- nce++;
- } else if (both) {
- double v1x = vx[2 * v1];
- double v2x = vx[2 * v2];
-
- if (fabs(v1x - v2x) > 0.5) {
- ce[nce] = i;
- nce++;
- }
- }
- }
-
- return nce;
-}
-
-graph_t *graph_create(lattice_t lattice, bound_t bound, uint_t L, bool dual, cholmod_common *c) {
- graph_t *g = (graph_t *)calloc(1, sizeof(graph_t));
- frame_t *f = frame_create(lattice, L, dual);
-
- g->L = L;
- g->boundary = bound;
- g->ne = f->ne;
-
- if (bound == TORUS_BOUND) {
- g->nv = f->nv;
- g->dnv = f->dnv;
- g->nb = 1;
-
- g->ev = f->ev;
- f->ev = NULL;
- g->dev = f->dev;
- f->dev = NULL;
- g->vx = f->vx;
- f->vx = NULL;
- g->dvx = f->dvx;
- f->dvx = NULL;
-
- // the boundary for the torus consists of a list of edges required to cut
- // the torus into a cylinder
- g->bi = (uint_t *)calloc(2, sizeof(uint_t));
- g->b = (uint_t *)malloc(g->ne * sizeof(uint_t));
- g->bi[1] = get_cut_edges(g->ne, g->ev, g->vx, false, g->b);
- g->bq = (bool *)calloc(g->ne, sizeof(bool));
- for (uint_t i = 0; i < g->bi[1]; i++) {
- g->bq[g->b[i]] = true;
- }
- } else {
- uint_t *ce = (uint_t *)malloc(g->ne * sizeof(uint_t));
- uint_t nce = 0;
-
- if (bound == CYLINDER_BOUND) {
- g->nb = 2;
- nce = get_cut_edges(f->ne, f->ev, f->vx, false, ce);
- } else {
- g->nb = 4;
- nce = get_cut_edges(f->ne, f->ev, f->vx, true, ce);
- }
-
- g->nv = f->nv;
- g->dnv = f->dnv;
- g->vx = (double *)malloc(2 * (f->nv + nce) * sizeof(double));
- g->dvx = (double *)malloc(2 * (f->dnv + nce) * sizeof(double));
- g->ev = f->ev;
- g->dev = f->dev;
- f->ev = NULL;
- f->dev = NULL;
- memcpy(g->vx, f->vx, 2 * f->nv * sizeof(double));
- memcpy(g->dvx, f->dvx, 2 * f->dnv * sizeof(double));
-
- uint_t nbv = 0;
- uint_t *bv = (uint_t *)calloc(f->nv, sizeof(uint_t));
- uint_t *dbv = (uint_t *)calloc(f->dnv, sizeof(uint_t));
- uint_t nside = 0;
- bool *side = (bool *)calloc(f->nv, sizeof(bool));
-
- for (uint_t i = 0; i < nce; i++) {
- uint_t cv1 = g->ev[2 * ce[i]];
- uint_t cv2 = g->ev[2 * ce[i] + 1];
- uint_t dv1 = g->dev[2 * ce[i]];
- uint_t dv2 = g->dev[2 * ce[i] + 1];
-
- uint_t cin = 1;
-
- if (bound == FREE_BOUND && (f->vx[2 * cv2] - f->vx[2 * cv1]) > 0.5) {
- cin = 0;
- }
-
- uint_t vin = f->vx[2 * cv1 + cin] < f->vx[2 * cv2 + cin] ? 0 : 1;
- uint_t dvin = f->dvx[2 * dv1 + cin] < f->dvx[2 * dv2 + cin] ? 0 : 1;
-
- if (bv[g->ev[2 * ce[i] + vin]] == 0) {
- nbv++;
- if (cin == 0) {
- nside++;
- side[g->ev[2 * ce[i] + vin]] = true;
- }
-
- bv[g->ev[2 * ce[i] + vin]] = g->nv;
-
- g->vx[2 * g->nv + cin] = 1 + f->vx[2 * g->ev[2 * ce[i] + vin] + cin];
- g->vx[2 * g->nv + !cin] = f->vx[2 * g->ev[2 * ce[i] + vin] + !cin];
- g->ev[2 * ce[i] + vin] = g->nv;
-
- g->nv++;
- } else {
- g->ev[2 * ce[i] + vin] = bv[g->ev[2 * ce[i] + vin]];
- }
- if (dbv[g->dev[2 * ce[i] + dvin]] == 0) {
- dbv[g->dev[2 * ce[i] + dvin]] = g->dnv;
-
- if (f->dvx[2 * g->dev[2 * ce[i] + dvin] + cin] < 0.5) {
- g->dvx[2 * g->dnv + cin] = 1 + f->dvx[2 * g->dev[2 * ce[i] + dvin] + cin];
- } else {
- g->dvx[2 * g->dnv + cin] = f->dvx[2 * g->dev[2 * ce[i] + dvin] + cin];
- }
- g->dvx[2 * g->dnv + !cin] = f->dvx[2 * g->dev[2 * ce[i] + dvin] + !cin];
- g->dev[2 * ce[i] + dvin] = g->dnv;
-
- g->dnv++;
- } else {
- g->dev[2 * ce[i] + dvin] = dbv[g->dev[2 * ce[i] + dvin]];
- }
- }
-
- g->bi = (uint_t *)calloc(1 + g->nb, sizeof(uint_t));
- g->b = (uint_t *)malloc(2 * nbv * sizeof(uint_t));
-
- g->bi[1] = ((int_t)nbv - (int_t)nside);
- g->bi[g->nb] = 2 * nbv;
-
- if (bound == FREE_BOUND) {
- g->bi[2] = 2 * ((int_t)nbv - (int_t)nside);
- g->bi[3] = 2 * (int_t)nbv - (int_t)nside;
- }
-
- uint_t n_ins0 = 0;
- uint_t n_ins1 = 0;
-
- g->nbi = (uint_t *)malloc(((int_t)g->nv - (int_t)g->bi[g->nb]) * sizeof(uint_t));
- g->bni = (uint_t *)malloc((g->nv + 1) * sizeof(uint_t));
- g->bq = (bool *)calloc(g->nv, sizeof(bool));
- uint_t n_tmp = 0;
-
- for (uint_t i = 0; i < f->nv; i++) {
- g->bni[i] = n_tmp;
- if (bv[i] != 0) {
- g->bq[i] = true;
- g->bq[bv[i]] = true;
- if (side[i]) {
- g->b[g->bi[2] + n_ins1] = i;
- g->b[g->bi[3] + n_ins1] = bv[i];
- n_ins1++;
- } else {
- g->b[g->bi[0] + n_ins0] = i;
- g->b[g->bi[1] + n_ins0] = bv[i];
- n_ins0++;
- }
- } else {
- g->nbi[n_tmp] = i;
- n_tmp++;
- }
- }
-
- for (uint_t i = 0; i < g->nv - f->nv + 1; i++) {
- g->bni[f->nv + i] = n_tmp;
- }
-
- free(bv);
- free(dbv);
- free(side);
- }
-
- g->vei = get_verts_to_edges_ind(g->nv, g->ne, g->ev);
- g->ve = get_verts_to_edges(g->nv, g->ne, g->ev, g->vei);
- g->dvei = get_verts_to_edges_ind(g->dnv, g->ne, g->dev);
- g->dve = get_verts_to_edges(g->dnv, g->ne, g->dev, g->dvei);
-
- g->voltcurmat = gen_voltcurmat(g->ne, g->nv, g->ev, c);
- uint_t num_spanning_edges;
- g->spanning_edges = get_spanning_edges(g->ne, g->ev, g->vx, 0.5, &num_spanning_edges);
- g->num_spanning_edges = num_spanning_edges;
-
- frame_free(f);
-
- return g;
-}
-
-
diff --git a/src/graph_free.c b/src/graph_free.c
deleted file mode 100644
index ee30c99..0000000
--- a/src/graph_free.c
+++ /dev/null
@@ -1,19 +0,0 @@
-
-#include "fracture.h"
-
-void graph_free(graph_t *g, cholmod_common *c) {
- free(g->ev);
- free(g->vei);
- free(g->ve);
- free(g->bi);
- free(g->b);
- free(g->vx);
- free(g->dvx);
- free(g->dev);
- free(g->dvei);
- free(g->dve);
- free(g->nbi);
- free(g->bq);
- CHOL_F(free_sparse)(&(g->voltcurmat), c);
- free(g);
-}
diff --git a/src/graph_genfunc.c b/src/graph_genfunc.c
deleted file mode 100644
index d396206..0000000
--- a/src/graph_genfunc.c
+++ /dev/null
@@ -1,19 +0,0 @@
-
-#include "fracture.h"
-
-double *genfunc_uniform(uint_t num, gsl_rng *r) {
- double *lattice = (double *)malloc(2 * num * sizeof(double));
-
- for (uint_t i = 0; i < 2 * num; i++) {
- lattice[i] = gsl_rng_uniform(r);
- }
-
- return lattice;
-}
-
-double *genfunc_hyperuniform(uint_t num, gsl_rng *r) {
- double *lattice = spheres(num, 0.8, 0.5, 0.9, 100, 100000);
-
- return lattice;
-}
-
diff --git a/src/net.c b/src/net.c
deleted file mode 100644
index a1047e2..0000000
--- a/src/net.c
+++ /dev/null
@@ -1,141 +0,0 @@
-
-#include "fracture.h"
-
-long double *get_thres(uint_t ne, double beta) {
- long double *thres = (long double *)malloc(ne * sizeof(long double));
- assert(thres != NULL);
-
- gsl_rng *r = gsl_rng_alloc(GSL_RAND_GEN);
- gsl_rng_set(r, rand_seed());
-
- for (uint_t i = 0; i < ne; i++) {
- thres[i] = rand_dist_pow(r, beta);
- }
-
- gsl_rng_free(r);
-
- return thres;
-}
-
-void net_notch(net_t *net, double notch_len, cholmod_common *c) {
- for (uint_t i = 0; i < net->graph->ne; i++) {
- uint_t v1, v2;
- double v1x, v1y, v2x, v2y, dy;
- bool crosses_center, not_wrapping, correct_length;
-
- v1 = net->graph->ev[2 * i];
- v2 = net->graph->ev[2 * i + 1];
-
- v1x = net->graph->vx[2 * v1];
- v1y = net->graph->vx[2 * v1 + 1];
- v2x = net->graph->vx[2 * v2];
- v2y = net->graph->vx[2 * v2 + 1];
-
- dy = v1y - v2y;
-
- crosses_center = (v1y >= 0.5 && v2y <= 0.5) || (v1y <= 0.5 && v2y >= 0.5);
- not_wrapping = fabs(dy) < 0.5;
- //correct_length = v1x + dx / dy * (v1y - 0.5) <= notch_len;
- correct_length = v1x < notch_len && v2x < notch_len;
-
- if (crosses_center && not_wrapping && correct_length) {
- break_edge(net, i, c);
- }
- }
-}
-
-net_t *net_create(const graph_t *g, double inf, double beta, double notch_len, bool vb, cholmod_common *c) {
- net_t *net = (net_t *)calloc(1, sizeof(net_t));
- assert(net != NULL);
-
- net->graph = g;
- net->num_broken = 0;
-
- net->fuses = (bool *)calloc(g->ne, sizeof(bool));
- assert(net->fuses != NULL);
-
- net->thres = get_thres(g->ne, beta);
- net->inf = inf;
-
- net->dim = g->nv;
-
- if (g->boundary == TORUS_BOUND) {
- net->nep = g->ne;
- net->evp = (uint_t *)malloc(2 * g->ne * sizeof(uint_t));
- memcpy(net->evp, g->ev, 2 * g->ne * sizeof(uint_t));
- } else {
- if (vb) {
- net->dim -= g->bi[g->nb];
- net->evp = (uint_t *)malloc(2 * g->ne * sizeof(uint_t));
- net->nep = 0;
- for (uint_t i = 0; i < g->ne; i++) {
- if (!(g->bq[g->ev[2*i]] || g->bq[g->ev[2*i+1]])) {
- net->evp[2*net->nep] = g->bni[g->ev[2*i]];
- net->evp[2*net->nep + 1] = g->bni[g->ev[2*i + 1]];
- net->nep++;
- }
- }
- } else {
- net->dim += 2;
- net->evp = (uint_t *)malloc(2 * (g->ne + g->bi[2]) * sizeof(uint_t));
- memcpy(net->evp, g->ev, 2 * g->ne * sizeof(uint_t));
- net->nep = g->ne + g->bi[2];
-
- for (uint_t i = 0; i < 2; i++) {
- for (uint_t j = 0; j < g->bi[i+1] - g->bi[i]; j++){
- net->evp[2 * (g->ne + g->bi[i] + j)] = g->b[g->bi[i] + j];
- net->evp[2 * (g->ne + g->bi[i] + j) + 1] = g->nv + i;
- }
- }
- }
- }
-
- net->voltage_bound = vb;
- net->boundary_cond = bound_set(g, vb, notch_len, c);
-
- net_notch(net, notch_len, c);
-
- {
- cholmod_sparse *laplacian = gen_laplacian(net, c);
- net->factor = CHOL_F(analyze)(laplacian, c);
- CHOL_F(factorize)(laplacian, net->factor, c);
- CHOL_F(free_sparse)(&laplacian, c);
- }
-
- return net;
-}
-
-net_t *net_copy(const net_t *net, cholmod_common *c) {
- net_t *net_copy = (net_t *)calloc(1, sizeof(net_t));
- assert(net_copy != NULL);
- memcpy(net_copy, net, sizeof(net_t));
-
- size_t fuses_size = (net->graph)->ne * sizeof(bool);
- net_copy->fuses = (bool *)malloc(fuses_size);
- assert(net_copy->fuses != NULL);
- memcpy(net_copy->fuses, net->fuses, fuses_size);
-
- size_t thres_size = (net->graph)->ne * sizeof(long double);
- net_copy->thres = (long double *)malloc(thres_size);
- assert(net_copy->thres != NULL);
- memcpy(net_copy->thres, net->thres, thres_size);
-
- size_t evp_size = 2 * net->nep * sizeof(uint_t);
- net_copy->evp = (uint_t *)malloc(thres_size);
- assert(net_copy->evp != NULL);
- memcpy(net_copy->evp, net->evp, evp_size);
-
- net_copy->boundary_cond = CHOL_F(copy_dense)(net->boundary_cond, c);
- net_copy->factor = CHOL_F(copy_factor)(net->factor, c);
-
- return net_copy;
-}
-
-void net_free(net_t *net, cholmod_common *c) {
- free(net->fuses);
- free(net->thres);
- CHOL_F(free_dense)(&(net->boundary_cond), c);
- CHOL_F(free_factor)(&(net->factor), c);
- free(net->evp);
- free(net);
-}
diff --git a/src/net_conductivity.c b/src/net_conductivity.c
deleted file mode 100644
index e9325bb..0000000
--- a/src/net_conductivity.c
+++ /dev/null
@@ -1,35 +0,0 @@
-
-#include "fracture.h"
-
-double net_conductivity(const net_t *net, const double *voltages) {
- if (net->voltage_bound) {
- // the voltage drop across the network is fixed to one with voltage
- // boundary conditions, so the conductivity is the total current flowing
- double tot_cur = 0;
- for (uint_t i = 0; i < net->graph->num_spanning_edges; i++) {
- uint_t e = net->graph->spanning_edges[i];
-
- if (!net->fuses[e]) {
- uint_t v1, v2, s1, s2;
- double v1y, v2y;
-
- v1 = net->graph->ev[2 * e];
- v2 = net->graph->ev[2 * e + 1];
-
- v1y = net->graph->vx[2 * v1 + 1];
- v2y = net->graph->vx[2 * v2 + 1];
-
- s1 = v1y < v2y ? v1 : v2;
- s2 = v1y < v2y ? v2 : v1;
-
- tot_cur += voltages[s1] - voltages[s2];
- }
- }
-
- return fabs(tot_cur);
- } else {
- // the current across the network is fixed to one with current boundary
- // conditions, so the conductivity is the inverse of the total voltage drop
- return 1 / fabs(voltages[net->graph->nv] - voltages[net->graph->nv + 1]);
- }
-}
diff --git a/src/net_currents.c b/src/net_currents.c
deleted file mode 100644
index 9bb2874..0000000
--- a/src/net_currents.c
+++ /dev/null
@@ -1,51 +0,0 @@
-
-#include "fracture.h"
-
-double *net_currents(const net_t *net, const double *voltages, cholmod_common *c) {
- uint_t ne = net->graph->ne;
- uint_t dim = net->graph->nv;
- cholmod_sparse *voltcurmat = net->graph->voltcurmat;
-
- cholmod_dense *x = CHOL_F(allocate_dense)(dim, 1, dim, CHOLMOD_REAL, c);
-
- double *tmp_x = x->x;
- x->x = (void *)voltages;
-
- cholmod_dense *y = CHOL_F(allocate_dense)(ne, 1, ne, CHOLMOD_REAL, c);
-
- double alpha[2] = {1, 0};
- double beta[2] = {0, 0};
- CHOL_F(sdmult)(voltcurmat, 0, alpha, beta, x, y, c);
-
- double *currents = (double *)malloc(ne * sizeof(double));
-
- for (int i = 0; i < ne; i++) {
- if (net->fuses[i]) {
- currents[i] = 0;
- } else {
- currents[i] = ((double *)y->x)[i];
- }
- }
-
- if (net->graph->boundary == TORUS_BOUND) {
- for (uint_t i = 0; i < net->graph->bi[1]; i++) {
- uint_t e = net->graph->b[i];
- uint_t v1 = net->graph->ev[2 * e];
- uint_t v2 = net->graph->ev[2 * e + 1];
- double v1y = net->graph->vx[2 * v1 + 1];
- double v2y = net->graph->vx[2 * v2 + 1];
-
- if (v1y > v2y) {
- currents[e] += 1;
- } else {
- currents[e] -= 1;
- }
- }
- }
-
- x->x = tmp_x;
- CHOL_F(free_dense)(&x, c);
- CHOL_F(free_dense)(&y, c);
-
- return currents;
-}
diff --git a/src/net_fracture.c b/src/net_fracture.c
deleted file mode 100644
index e7f18fc..0000000
--- a/src/net_fracture.c
+++ /dev/null
@@ -1,67 +0,0 @@
-
-#include "fracture.h"
-
-uint_t get_next_broken(net_t *net, double *currents, double cutoff) {
- uint_t max_pos = UINT_MAX;
- long double max_val = 0;
-
- for (uint_t i = 0; i < net->graph->ne; i++) {
- long double current = fabs(currents[i]);
- bool broken = net->fuses[i];
-
- if (!broken && current > cutoff) {
- long double val = current / net->thres[i];
-
- if (val > max_val) {
- max_val = val;
- max_pos = i;
- }
- }
- }
-
- if (max_pos == UINT_MAX) {
- printf("GET_NEXT_BROKEN: currents is zero or NaN, no max_val found\n");
- exit(EXIT_FAILURE);
- }
-
- return max_pos;
-}
-
-data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff) {
- data_t *data = data_create(net->graph->ne);
-
- uint_t n = 0;
- while (true) {
- n++;
- double *voltages = net_voltages(net, c);
- double *currents = net_currents(net, voltages, c);
-
- double conductivity = net_conductivity(net, voltages);
-
- if (conductivity < cutoff) {
- free(voltages);
- free(currents);
- break;
- }
-
- uint_t last_broke = get_next_broken(net, currents, cutoff);
-
- long double sim_current;
-
- if (net->voltage_bound) {
- sim_current = conductivity;
- } else {
- sim_current = 1;
- }
-
- data_update(data, last_broke, fabsl(sim_current * (net->thres)[last_broke] / ((long double)currents[last_broke])), conductivity);
-
- free(voltages);
- free(currents);
-
- break_edge(net, last_broke, c);
- }
-
- return data;
-}
-
diff --git a/src/net_voltages.c b/src/net_voltages.c
deleted file mode 100644
index c3537a5..0000000
--- a/src/net_voltages.c
+++ /dev/null
@@ -1,40 +0,0 @@
-
-#include "fracture.h"
-
-double *net_voltages(const net_t *net, cholmod_common *c) {
- cholmod_dense *b = net->boundary_cond;
- cholmod_factor *factor = net->factor;
-
- cholmod_dense *x = CHOL_F(solve)(CHOLMOD_A, factor, b, c);
-
- if (((double *)x->x)[0] != ((double *)x->x)[0]) {
- printf("GET_VOLTAGE: value is NaN\n");
- exit(EXIT_FAILURE);
- }
-
- double *t_voltages = (double *)x->x;
- x->x = NULL;
- CHOL_F(free_dense)(&x, c);
-
- graph_t *g = net->graph;
-
- if (g->boundary == TORUS_BOUND) {
- return t_voltages;
- } else if (net->voltage_bound) {
- double *voltages = (double *)malloc(g->nv * sizeof(double));
- for (uint_t i = 0; i < g->nv - g->bi[g->nb]; i++) {
- voltages[g->nbi[i]] = t_voltages[i];
- }
- for (uint_t i = 0; i < 2; i++) {
- for (uint_t j = 0; j < g->bi[i + 1] - g->bi[i]; j++) {
- voltages[g->b[g->bi[i] + j]] = 1 - i;
- }
- }
-
- free(t_voltages);
- return voltages;
- } else {
- return t_voltages;
- }
-}
-
diff --git a/src/rand.c b/src/rand.c
deleted file mode 100644
index 75722ac..0000000
--- a/src/rand.c
+++ /dev/null
@@ -1,23 +0,0 @@
-
-#include "fracture.h"
-
-unsigned long int rand_seed() {
- FILE *f = fopen("/dev/urandom", "r");
- unsigned long int seed;
- fread(&seed, sizeof(unsigned long int), 1, f);
- fclose(f);
- return seed;
-}
-
-long double rand_dist_pow(const gsl_rng *r, double beta) {
- long double x = 0;
-
- // underflow means that for very small beta x is sometimes identically zero,
- // which causes problems
- while (x == 0.0) {
- long double y = logl(gsl_rng_uniform_pos(r)) / beta;
- x = expl(y);
- }
-
- return x;
-}
diff --git a/src/spheres_poly/box.cpp b/src/spheres_poly/box.cpp
deleted file mode 100644
index 716090d..0000000
--- a/src/spheres_poly/box.cpp
+++ /dev/null
@@ -1,1157 +0,0 @@
-/*
- 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 <iostream>
-#include <math.h>
-#include <stdlib.h>
-#include <time.h>
-#include <iomanip>
-
-//==============================================================
-//==============================================================
-// 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<DIM>[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<N; i++) // initialize binlist to -1
- binlist[i] = -1;
-
- time(&start);
-}
-
-
-//==============================================================
-// Destructor
-//==============================================================
-box::~box()
-{
- delete[] s;
- delete[] binlist;
- delete[] x;
-}
-
-
-//==============================================================
-// ReadFile
-//==============================================================
-void box::ReadPositions(const char* filename)
-{
- // open file to read in arrays
- std::ifstream infile(filename);
-
- infile.ignore(256, '\n'); // ignore the dim line
- infile.ignore(256, '\n'); // ignore the #sphere 1 line
- infile.ignore(256, '\n'); // ignore the #sphere line
- infile.ignore(256, '\n'); // ignore the diameter line
- infile.ignore(1000, '\n'); // ignore the 100 010 001 line
- infile.ignore(256, '\n'); // ignore the T T T line
-
- for (int i=0; i<N; i++)
- {
- infile >> 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<DIM; k++)
- infile >> 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<N; i++)
- {
- CreateSphere(Ncurrent);
- Ncurrent++;
- }
- if (Ncurrent != N)
- std::cout << "problem! only made " << Ncurrent << " out of " << N << " desired spheres" << std::endl;
-
- VelocityGiver(temp);
- SetInitialEvents();
-}
-
-
-//==============================================================
-// Creates a sphere of random radius r at a random unoccupied position
-//==============================================================
-void box::CreateSphere(int Ncurrent)
-{
- int keeper; // boolean variable: 1 means ok, 0 means sphere already there
- int counter = 0; // counts how many times sphere already exists
- vector<DIM> 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<DIM; k++)
- xrand[k] = ((double)random()/(double)RAND_MAX)*SIZE;
-
- for (int i=0; i<Ncurrent; i++) // check if overlapping other spheres
- {
- d=0.;
- if (hardwallBC)
- {
- for (int k=0; k<DIM; k++)
- {
- d += (xrand[k] - s[i].x[k])*(xrand[k] - s[i].x[k]);
- }
- }
- else
- {
- for (int k=0; k<DIM; k++)
- {
- if ((xrand[k] - s[i].x[k])*(xrand[k] - s[i].x[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<DIM; k++) // check if overlapping wall
- {
- if ((xrand[k] <= radius) || (SIZE - xrand[k] <= radius)) // touching wall
- {
- keeper = 0;
- counter++;
- break;
- }
- }
- }
-
- if (keeper == 1)
- break;
-
- if (counter >= 1000)
- {
- std::cout << "counter >= 1000" << std::endl;
- exit(-1);
- }
- }
-
- // now convert xrand into index vector for cells
- vector<DIM,int> cell;
- cell = vector<DIM>::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<N; i++)
- {
- // now convert x into index vector for cells
- vector<DIM,int> cell;
- cell = vector<DIM>::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; i<N; i++)
- {
- for (int k=0; k<DIM; k++)
- {
- if (T==0.)
- s[i].v[k] = 0.;
- else
- s[i].v[k] = Velocity(T);
- }
- }
-}
-
-
-//==============================================================
-// Velocity, gives a single velocity from Max/Boltz dist.
-//==============================================================
-double box::Velocity(double T)
-{
- double rand; // random number between -0.5 and 0.5
- double sigmasquared = T; // Assumes M = mass of sphere = 1
- double sigma = sqrt(sigmasquared); // variance of Gaussian
- double stepsize = 1000.; // stepsize for discretization of integral
- double vel = 0.0; // velocity
- double dv=sigma/stepsize;
- double p=0.0;
-
- rand = (double)random() / (double)RAND_MAX - 0.5;
- if(rand < 0)
- {
- rand = -rand;
- dv = -dv;
- }
-
- while(fabs(p) < rand) // integrate until the integral equals rand
- {
- p += dv * 0.39894228 * exp(-vel*vel/(2.*sigmasquared))/sigma;
- vel += dv;
- }
- return vel;
-}
-
-
-//==============================================================
-// Finds next events for all spheres..do this once at beginning
-//==============================================================
-void box::SetInitialEvents()
-{
- for (int i=0; i<N; i++) // set all events to checks
- {
- event e(gtime, i, INF);
- s[i].nextevent = e;
- h.insert(i);
- }
-}
-
-
-//==============================================================
-// Finds next event for sphere i
-//==============================================================
-event box::FindNextEvent(int i)
-{
- double outgrowtime;
- outgrowtime = (.5/ngrids - (s[i].r+s[i].gr*gtime))/s[i].gr + gtime;
-
- event t = FindNextTransfer(i);
- event c = FindNextCollision(i);
-
- if ((outgrowtime < c.time)&&(outgrowtime < t.time)&&(ngrids>1))
- {
- 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<DIM> xi = s[i].x + s[i].v*(gtime - s[i].lutime);
- vector<DIM> vi = s[i].v;
-
- for (int k=0; k<DIM; k++)
- {
- double newtime;
- if (vi[k]==0.)
- newtime= dblINF;
- else if (vi[k]>0) // 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<ttime)
- {
- wallindex = k+1;
- ttime = newtime;
- }
- }
- else if (vi[k]<0) // will hit left wall
- {
- if ((hardwallBC)&&(s[i].cell[k] == 0))
- newtime = ((double)(s[i].cell[k])*SIZE/((double)(ngrids))
- - (xi[k]-(s[i].r+s[i].gr*gtime)))/(vi[k]-s[i].gr);
- else
- newtime = ((double)(s[i].cell[k])*SIZE/((double)(ngrids))
- - xi[k])/(vi[k]);
-
- if (newtime<ttime)
- {
- wallindex = -(k+1);
- ttime = newtime;
- }
- }
- }
-
-
- if (ttime < 0)
- ttime = 0;
- // make the event and return it
- event e = event(ttime+gtime,i,wallindex+DIM+N+1);
- return e;
-}
-
-
-//==============================================================
-// Check all nearest neighbor cells for collision partners
-//==============================================================
-void box::ForAllNeighbors(int i, vector<DIM,int> vl, vector<DIM,int> vr,
- neighbor& operation)
-{
- vector<DIM,int> cell = s[i].cell;
-
- // now iterate through nearest neighbors
- vector<DIM, int> offset; // nonnegative neighbor offset
- vector<DIM, int> pboffset; // nearest image offset
-
- vector<DIM,int> grid;
-
- int ii ;
-
- grid=vl;
- while(1)
- {
- //if (vr[0] > 1)
- //std::cout << grid << "..." << cell+grid << "\n";
- for(int k=0; k<DIM; k++)
- {
- offset[k]=grid[k]+ngrids; // do this so no negatives
- if (cell[k]+grid[k]<0) //out of bounds to left
- pboffset[k] = -1;
- else if (cell[k]+grid[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;ii++)
- {
- grid[ii] += 1;
- if(grid[ii]<=vr[ii]) break;
- grid[ii]=vl[ii];
- }
- if(ii>=DIM) break;
- }
-}
-
-
-//==============================================================
-// PredictCollision
-//==============================================================
-void box::PredictCollision(int i, int j, vector<DIM, int> pboffset,
- double& ctime, int& cpartner,
- vector<DIM, int>& 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<DIM, int> vl, vr;
-
- for (int k=0; k<DIM; k++) // check all nearest neighbors
- {
- vl[k] = -1;
- vr[k] = 1;
- }
-
- ForAllNeighbors(i,vl,vr,cc);
-
- event e;
- if (cc.cpartner == i) // found no collisions in neighboring cells
- {
- if (cc.ctime != dblINF)
- std::cout << "ctime != dblINF" << std::endl;
- e = event(dblINF,i,INF); // give check at double INF
- }
- else
- e = event(cc.ctime,i,cc.cpartner,cc.cpartnerpboffset);
-
- return e;
-}
-
-
-//==============================================================
-// Calculates collision time between i and image of j using quadratic formula
-//==============================================================
-double box::CalculateCollision(int i, int j, vector<DIM> pboffset)
-{
- if ((hardwallBC)&&(pboffset.norm_squared() > 1E-12))
- {
- return dblINF;
- }
- else
- {
- // calculate updated position and velocity of i and j
- vector<DIM> xi = s[i].x + s[i].v*(gtime - s[i].lutime);
- vector<DIM> vi = s[i].v;
- vector<DIM> xj = s[j].x + pboffset*SIZE + s[j].v*(gtime - s[j].lutime);
- vector<DIM> 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<DIM>::norm_squared(vi - vj) - (s[i].gr+s[j].gr)*(s[i].gr+s[j].gr);
- B = vector<DIM>::dot(xi - xj, vi - vj) - (r_now_i+r_now_j)*(s[i].gr+s[j].gr);
- C = vector<DIM>::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<N)) // collision!
- {
- ncollisions++;
- //std::cout << "collision between " << e.i << " and " << e.j << " at time " << e.time << std::endl;
- Collision(e);
- f = FindNextEvent(i);
- s[i].nextevent = f;
- h.downheap(1);
- if (f.time < e.time)
- {
- std::cout << "error, replacing event with < time" << std::endl;
- exit(-1);
- }
-
- // make sure collision was symmetric and give j a check
- if ((s[e.j].nextevent.j != i)||(s[e.j].nextevent.time != gtime))
- {
- std::cout << "error collisions not symmetric" << std::endl;
- std::cout << "collision between " << e.i << " and " << e.j << " at time " << e.time << std::endl;
- std::cout << "but " << e.j << " thinks it has " << s[e.j].nextevent.j<< " " << s[e.j].nextevent.time << std::endl;
- exit(-1);
- }
- else // give j a check
- s[e.j].nextevent.j = INF;
- }
- else if (e.j==INF) // check!
- {
- nchecks++;
- //std::cout << "check for " << e.i << " at time " << e.time << std::endl;
- f = FindNextEvent(i);
- s[i].nextevent = f;
- h.downheap(1);
- }
- else if (e.j==INF-1) // sphere outgrowing unit cell, decrease ngrids!
- {
- gtime = e.time;
- Synchronize(false);
- ngrids = ngrids - 1;
-// std::cout << "need to reduce ngrids to " << ngrids << std::endl;
- ChangeNgrids(ngrids);
- h.downheap(1);
- }
- else // transfer!
- {
- ntransfers++;
- //std::cout << "transfer for " << e.i << " at time " << e.time << std::endl;
- Transfer(e);
- f = FindNextEvent(i);
- s[i].nextevent = f;
- h.downheap(1);
- //r = FindNextEvent(i, e.j-N-DIM-1);
- //if (f.time <= e.time)
- if (f.time < e.time)
- {
- std::cout << "error after transfer, replacing new event with < time" << " " << std::endl;
- std::cout << std::setprecision(16) << "e.time= " << e.time << ", f.time= " << f.time << ", f.i= " << f.i << ", f.j= " << f.j << "e.i= " << e.i << ", e.j= " << e.j << std::endl;
- std::cout << std::setprecision(16) << "difference= " << e.time - f.time << std::endl;
- exit(-1);
- }
- }
-}
-
-
-//==============================================================
-// Processes a collision
-//=============================================================
-void box::Collision(event e)
-{
- double ctime = e.time;
- int i = e.i;
- int j = e.j;
- vector<DIM,int> v = e.v; // virtual image
- gtime = ctime;
-
- // Update positions and cells of i and j to ctime
- s[i].x += s[i].v*(gtime-s[i].lutime);
- s[j].x += s[j].v*(gtime-s[j].lutime);
-
- // Check to see if a diameter apart
- double r_sum = s[i].r + s[j].r + (s[i].gr+s[j].gr)*gtime;
- double distance = vector<DIM>::norm_squared(s[i].x - s[j].x- v.Double()*SIZE) - r_sum*r_sum;
- if (distance*distance > 10.*DBL_EPSILON)
- std::cout << "overlap " << distance << std::endl;
-
- s[i].lutime = gtime;
- s[j].lutime = gtime;
-
- vector<DIM,double> vipar; // parallel comp. vi
- vector<DIM,double> vjpar; // parallel comp. vj
- vector<DIM,double> viperp; // perpendicular comp. vi
- vector<DIM,double> vjperp; // perpendicular comp. vj
-
- // make unit vector out of displacement vector
- vector<DIM,double> 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<DIM>::dot(s[i].v, dhat);
- vjpar = dhat*vector<DIM>::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<DIM>::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<DIM,int> 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 (j<N+DIM+1) // left wall
- {
- k = -j+N+DIM;
- celli[k] = s[i].cell[k] - 1;
-
- if (hardwallBC)
- {
- // if in left-most cell, reflect back
- if (s[i].cell[k] == 0)
- s[i].v[k] = -s[i].v[k] + 4.*s[i].gr;
- else
- {
- if (ngrids>1)
- 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<DIM,int>& 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<N; i++)
- std::cout << i << " " << s[i].x << " " << s[i].v << " " << s[i].cell << std::endl;
-}
-
-
-//==============================================================
-// Update positions...purely for graphical display
-//==============================================================
-void box::TrackPositions()
-{
- for (int i=0; i<N; i++)
- x[i] = s[i].x + s[i].v*(gtime-s[i].lutime);
-}
-
-
-//==============================================================
-// Computes the total energy
-//==============================================================
-double box::Energy()
-{
- double E=0;
- for (int i=0; i<N; i++)
- E += 0.5*s[i].m*s[i].v.norm_squared();
-
- return E/N;
-}
-
-
-//==============================================================
-// Calculates the packing fraction
-//==============================================================
-double box::PackingFraction()
-{
- double rfactor = 0.;
- for (int i=0; i<N; i++)
- {
- rfactor += pow(s[i].r + gtime*s[i].gr, DIM);
- }
- double v = (rfactor*pow(sqrt(PI), DIM))/(exp(lgamma(1.+((double)(DIM))/2.)));
- return v/(pow(SIZE, DIM));
-}
-
-
-//==============================================================
-// Checks to make sure all sphere diameters are less than dimension
-// of unit cell
-//==============================================================
-int box::CheckSphereDiameters()
-{
- int offender = 0;
- for (int i=0; i<N; i++)
- {
- if (s[i].r*2 > 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<N; i++) // initialize binlist to -1
- binlist[i] = -1;
- AssignCells();
- for (int i=0; i<N; i++)
- s[i].nextevent = event(0., i, INF);
- Process(N);
-}
-
-
-//==============================================================
-// Calculates the optimal ngrids for the initial configuration
-// and assumes that ngrids gets updated (reduced) as the packing
-// proceeds
-//==============================================================
-int box::Optimalngrids2(double currentradius)
-{
- return (int)(1./(2.*currentradius));
-}
-
-
-//==============================================================
-// Calculates the optimal ngrids, assuming ngrids is not updated
-// automatically and is very conservative
-//==============================================================
-int box::Optimalngrids()
-{
- double maxr;
-
- maxr = pow(exp(lgamma(1.+((double)(DIM))/2.))*maxpf/
- (N*(bidispersityfraction + (1.-bidispersityfraction)*
- pow(bidispersityratio, DIM))),
- 1./DIM)/sqrt(PI);
-
- return (int)(1./(2.*maxr));
-}
-
-
-//==============================================================
-// Processes n events
-//==============================================================
-void box::Process(int n)
-{
- double deltat = gtime;
- for (int i=0; i<n; i++)
- {
- ProcessEvent();
- }
- pf = PackingFraction(); // packing fraction
- deltat = gtime - deltat;
- double oldenergy = energy;
- energy = Energy(); // kinetic energy
-
- energychange = ((oldenergy - energy)/oldenergy)*100; // percent change in energy
-
- if (deltat != 0.)
- {
- pressure = 1+xmomentum/(2.*energy*N*deltat);
- collisionrate = ((double)(ncollisions))/deltat;
- }
-
- // reset to 0
- ncollisions = 0;
- ntransfers = 0;
- nchecks = 0;
- xmomentum = 0.;
- ncycles++;
-}
-
-
-//==============================================================
-// Prints statistics for n events
-//==============================================================
-void box::PrintStatistics()
-{
- std::cout << "packing fraction = " << pf << std::endl;
- std::cout << "gtime = " << gtime << std::endl;
- std::cout << "total time = " << rtime+gtime << std::endl;
- std::cout << "kinetic energy = " << energy << std::endl;
- std::cout << "total # events = " << neventstot << std::endl;
- std::cout << "# events = " << ncollisions+ntransfers+nchecks << ", # collisions = " << ncollisions << ", # transfers = " << ntransfers << ", # checks =" << nchecks << std::endl;
- std::cout << "growthrate = " << growthrate << std::endl;
- std::cout << "collisionrate = " << collisionrate << std::endl;
- std::cout << "reduced pressure = " << pressure << std::endl;
- std::cout << "-----------------" << std::endl;
-}
-
-
-//==============================================================
-// Updates spheres to gtime, synchronizes, and can change growth rate
-//==============================================================
-void box::Synchronize(bool rescale)
-{
- double vavg = sqrt(2.*M*energy);
-
- for (int i=0; i<N; i++)
- {
- s[i].x = s[i].x + s[i].v*(gtime-s[i].lutime);
- s[i].nextevent.time -= gtime;
-
- if (s[i].nextevent.time < 0.)
- std::cout << "error, event times negative after synchronization" << std::endl;
- if (rescale == true) // give everyone checks
- {
- s[i].nextevent = event(0., i, INF);
- s[i].v /= vavg;
- }
-
- s[i].lutime = 0.;
- s[i].r += gtime*s[i].gr;
- }
-
- //r += gtime*growthrate; // r defined at gtime = 0
- rtime += gtime;
- gtime = 0.;
-
- if (rescale == true)
- Process(N);
-}
-
-
-//==============================================================
-// Run time
-//==============================================================
-void box::RunTime()
-{
- time(&end);
- std::cout << "run time = " << difftime(end, start) << std::endl;
-}
-
-
-//==============================================================
-// Write configuration
-//==============================================================
-void box::WriteConfiguration(double *data)
-{
- int count1;
-
- if (gtime != 0.) // synchronize spheres if not currently synchronized
- Synchronize(false);
-
-
- for (int i=0; i<N; i++) // output radius and positions
- {
- for (int k=0; k<DIM; k++)
- data[DIM * i + k] = s[i].x[k];
- }
-
-}
diff --git a/src/spheres_poly/box.h b/src/spheres_poly/box.h
deleted file mode 100644
index 69418f4..0000000
--- a/src/spheres_poly/box.h
+++ /dev/null
@@ -1,169 +0,0 @@
-/*
- 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.
-*/
-
-//-----------------------------------------------------------------------------
-// Box maker
-//---------------------------------------------------------------------------
-
-#ifndef BOX_H
-#define BOX_H
-
-
-#include <vector>
-#include <math.h>
-
-#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<DIM, int>& 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<DIM, int>, vector<DIM,int>, neighbor&);
- void PredictCollision(int i, int j, vector<DIM, int> pboffset,
- double& ctime, int& cpartner,
- vector<DIM, int>& cpartnerpboffset);
- double CalculateCollision(int i, int j, vector<DIM> 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<DIM,int>& 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<DIM, int> cells; // array that keeps track of spheres in each cell
- int *binlist; // linked-list for cells array
- heap h; // event heap
- vector<DIM> *x; // positions of spheres.used for graphics
-};
-
-
-//---------------------------------------------------------------------------
-// Predicts collisions, inherits neighbor operation
-//---------------------------------------------------------------------------
-class collision : public neighbor
-{
- public:
-
- box *b;
- double ctime;
- int cpartner;
- vector<DIM,int> cpartnerpboffset;
-
- public:
- collision(int i_i, box *b);
-
- virtual void Operation(int j, vector<DIM, int>& pboffset);
-};
-
-#endif
diff --git a/src/spheres_poly/event.cpp b/src/spheres_poly/event.cpp
deleted file mode 100644
index c8c104d..0000000
--- a/src/spheres_poly/event.cpp
+++ /dev/null
@@ -1,65 +0,0 @@
-#include "event.h"
-
-//==============================================================
-//==============================================================
-// Class Event
-//==============================================================
-//==============================================================
-
-
-//==============================================================
-// Constructor
-//==============================================================
-event::event(double time_i, int i_i, int j_i, vector<DIM,int> 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
deleted file mode 100644
index 55dd1fc..0000000
--- a/src/spheres_poly/event.h
+++ /dev/null
@@ -1,58 +0,0 @@
-#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<DIM,int> 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<DIM,int> 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
deleted file mode 100644
index 10430ed..0000000
--- a/src/spheres_poly/grid_field.h
+++ /dev/null
@@ -1,148 +0,0 @@
-/*
- 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<int D, class T>
-class grid_field {
-
- public:
- int elements;
-
- private:
- T* f;
- vector<D, int> size; // number of grid points for each dimension
- vector<D, int> offset;
-
- public:
-
- grid_field();
- grid_field(const vector<D, int>&);
- ~grid_field();
-
- T& get(const vector<D, int>&);
-
- vector<D, int> get_size() const;
- void set_size(const vector<D, int>&);
- void set_size(const int);
- void initialize(const int i);
-};
-
-
-// grid_field
-// ~~~~~~~~~~~~
-template<int D, class T>
-grid_field<D, T>::grid_field()
- : f(0), elements(0)
-{
-}
-
-
-// grid_field
-// ~~~~~~~~~~~~
-template<int D, class T>
-grid_field<D, T>::grid_field(const vector<D, int>& s)
- : f(0)
-{
- set_size(s);
-}
-
-
-// ~grid_field
-// ~~~~~~~~~~~~~
-template <int D, class T>
-grid_field<D, T>::~grid_field()
-{
- if(f != 0)
- delete[] f;
-}
-
-
-// get_size
-// ~~~~~~~~
-template<int D, class T>
-inline vector<D, int> grid_field<D, T>::get_size() const
-{
- return size;
-}
-
-
-// set_size
-// ~~~~~~~~
-template<int D, class T>
-void grid_field<D, T>::set_size(const vector<D, int>& s)
-{
- if(f != 0)
- delete[] f;
-
- size = s;
-
- elements = 1;
- for(int i=0; i<D; i++) {
- offset[i] = elements;
- elements *= size.x[i];
- }
-
- f = new T[elements];
-}
-
-
-// set_size
-// ~~~~~~~~
-template<int D, class T>
-void grid_field<D, T>::set_size(const int s)
-{
- vector<D, int> square;
-
- for(int k=0; k<D; k++)
- square[k] = s;
-
- set_size(square);
-}
-
-
-// get
-// ~~~
-template<int D, class T>
-inline T& grid_field<D, T>::get(const vector<D, int>& pos)
-{
- int p=0;
- for(int i=0; i<D; i++)
- p += pos.x[i]*offset[i];
-
- return f[p];
-}
-
-
-// initialize
-// ~~~
-template<int D, class T>
-void grid_field<D, T>::initialize(const int value)
-{
- for(int i=0; i<elements; i++)
- f[i] = value;
-}
-
-
-
-#endif
diff --git a/src/spheres_poly/heap.cpp b/src/spheres_poly/heap.cpp
deleted file mode 100644
index c83865c..0000000
--- a/src/spheres_poly/heap.cpp
+++ /dev/null
@@ -1,186 +0,0 @@
-#include "heap.h"
-#include "event.h"
-#include <iostream>
-#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
deleted file mode 100644
index 91180e3..0000000
--- a/src/spheres_poly/heap.h
+++ /dev/null
@@ -1,42 +0,0 @@
-//---------------------------------------------------------------------------
-// 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
deleted file mode 100644
index 920c099..0000000
--- a/src/spheres_poly/neighbor.cpp
+++ /dev/null
@@ -1,40 +0,0 @@
-#include "vector.h"
-#include <iostream>
-#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<DIM, int>& pboffset)
-{
- b->PredictCollision(i, j, pboffset, ctime, cpartner, cpartnerpboffset);
-}
-
diff --git a/src/spheres_poly/sphere.cpp b/src/spheres_poly/sphere.cpp
deleted file mode 100644
index 0fbbd68..0000000
--- a/src/spheres_poly/sphere.cpp
+++ /dev/null
@@ -1,66 +0,0 @@
-#include "box.h"
-#include <iostream>
-#include <math.h>
-#include <stdlib.h>
-#include <time.h>
-
-#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<DIM> x_i, vector<DIM, int> 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
deleted file mode 100644
index 28c64b4..0000000
--- a/src/spheres_poly/sphere.h
+++ /dev/null
@@ -1,44 +0,0 @@
-#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<DIM> x, vector<DIM, int> 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<DIM, int> cell; // cell that it belongs to
- vector<DIM, double> x; // position
- vector<DIM, double> 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
deleted file mode 100644
index 9b14b01..0000000
--- a/src/spheres_poly/spheres.cpp
+++ /dev/null
@@ -1,64 +0,0 @@
-//===========================================================
-//===========================================================
-//===========================================================
-//
-// Molecular dynamics simulation of hardspheres
-//
-//===========================================================
-//===========================================================
-//===========================================================
-
-#include <iostream>
-#include <math.h>
-#include <fstream>
-#include <vector>
-#include <time.h>
-#include <string.h>
-#include <stdlib.h>
-#include <cstdlib>
-
-#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
deleted file mode 100644
index 44d8052..0000000
--- a/src/spheres_poly/spheres.h
+++ /dev/null
@@ -1,4 +0,0 @@
-
-#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
deleted file mode 100644
index 9ee481f..0000000
--- a/src/spheres_poly/vector.h
+++ /dev/null
@@ -1,471 +0,0 @@
-#ifndef VECTOR_H
-#define VECTOR_H
-
-#include <iostream>
-#include <fstream>
-
-#define DIM 2
-
-template <int D, typename T=double>
-class vector {
-
- public:
- T x[D];
-
- public:
-
- vector();
- vector(const T[D]);
- vector(const vector&);
- ~vector();
-
- vector<D, T>& operator+=(const vector<D, T>&);
- vector<D, T>& operator-=(const vector<D, T>&);
- vector<D, T>& operator*=(const T);
- vector<D, T>& operator/=(const T);
- vector<D, T> operator+(const vector<D, T>&) const;
- vector<D, T> operator-(const vector<D, T>&) const;
- vector<D, T> operator*(const T) const;
- vector<D, T> operator/(const T) const;
- vector<D, T> operator%(const T) const;
- bool operator==(const vector<D, T> &a) const;
-
- vector<D, int> integer() const;
- vector<D, double> Double() const;
- static vector<D, int> integer(const vector<D, T>&);
- static vector<D, double> Double(const vector<D, T>&);
-
- T& operator[](const unsigned int);
-
- double dot(const vector<D, T>&) const;
- static double dot(const vector<D, T>&, const vector<D, T>&);
-
- double norm_squared() const;
- static double norm_squared(const vector<D, T>&);
-
- void read(std::ifstream&);
- void write(std::ofstream&) const;
-};
-
-
-template <int D, typename T>
-std::ostream& operator<<(std::ostream&, const vector<D, T>&);
-
-
-// constructor
-// ~~~~~~~~~~~
-template <int D, typename T>
-vector<D, T>::vector()
-{
- for(int k=0; k<D; k++)
- x[k] = 0;
-}
-
-template <int D, typename T>
-vector<D, T>::vector(const T x_i[D])
-{
- for(int k=0; k<D; k++)
- x[k] = x_i[k];
-}
-
-template <int D, typename T>
-vector<D, T>::vector(const vector<D, T> &v)
-{
- for(int k=0; k<D; k++)
- x[k] = v.x[k];
-}
-
-
-// destructor
-// ~~~~~~~~~~
-template <int D, typename T>
-vector<D, T>::~vector()
-{
-}
-
-
-// +=
-// ~~
-template <int D, typename T>
-inline vector<D, T>& vector<D, T>::operator+=(const vector<D, T> &v)
-{
- for(int k=0; k<D; k++)
- x[k] += v.x[k];
-
- return *this;
-}
-
-
-// -=
-// ~~
-template <int D, typename T>
-inline vector<D, T>& vector<D, T>::operator-=(const vector<D, T> &v)
-{
- for(int k=0; k<D; k++)
- x[k] -= v.x[k];
-
- return *this;
-}
-
-
-// *=
-// ~~
-template <int D, typename T>
-inline vector<D, T>& vector<D, T>::operator*=(const T s)
-{
- for(int k=0; k<D; k++)
- x[k] *= s;
-
- return *this;
-}
-
-// /=
-// ~~
-template <int D, typename T>
-inline vector<D, T>& vector<D, T>::operator/=(const T s)
-{
- for(int k=0; k<D; k++)
- x[k] /= s;
-
- return *this;
-}
-
-
-// +
-// ~
-template <int D, typename T>
-inline vector<D, T> vector<D, T>::operator+(const vector<D, T> &a) const
-{
- vector<D, T> c;
-
- for(int k=0; k<D; k++)
- c.x[k] = x[k] + a.x[k];
-
- return c;
-}
-
-
-// -
-// ~
-template <int D, typename T>
-inline vector<D, T> vector<D, T>::operator-(const vector<D, T> &a) const
-{
- vector<D, T> c;
-
- for(int k=0; k<D; k++)
- c[k] = x[k] - a.x[k];
-
- return c;
-}
-
-
-// *
-// ~
-template <int D, typename T>
-inline vector<D, T> vector<D, T>::operator*(const T s) const
-{
- vector<D, T> c;
-
- for(int k=0; k<D; k++)
- c[k] = x[k] * s;
-
- return c;
-}
-
-
-// /
-// ~
-template <int D, typename T>
-inline vector<D, T> vector<D, T>::operator/(const T s) const
-{
- vector<D, T> c;
-
- for(int k=0; k<D; k++)
- c[k] = x[k] / s;
-
- return c;
-}
-
-
-// ==
-// ~
-template <int D, typename T>
-inline bool vector<D, T>::operator==(const vector<D, T> &a) const
-{
- for(int k=0; k<D; k++)
- {
- if (!(x[k]==a.x[k]))
- return false;
- }
- return true;
-}
-
-
-// %
-// ~
-template <int D, typename T>
-inline vector<D, T> vector<D, T>::operator%(const T s) const
-{
- vector<D, T> c;
-
- for(int k=0; k<D; k++)
- c[k] = x[k] % s;
-
- return c;
-}
-
-
-// integer
-// ~~~~~~~
-template <int D, typename T>
-inline vector<D, int> vector<D, T>::integer() const
-{
- vector<D, int> c;
-
- for(int k=0; k<D; k++)
- c[k] = (int)x[k];
-
- return c;
-}
-
-template <int D, typename T>
-inline vector<D, int> vector<D, T>::integer(const vector<D, T>& v)
-{
- return v.integer();
-}
-
-
-// double
-// ~~~~~~~
-template <int D, typename T>
-inline vector<D, double> vector<D, T>::Double() const
-{
- vector<D, double> c;
-
- for(int k=0; k<D; k++)
- c[k] = (double)x[k];
-
- return c;
-}
-
-template <int D, typename T>
-inline vector<D, double> vector<D, T>::Double(const vector<D, T>& v)
-{
- return v.Double();
-}
-
-
-
-// []
-// ~~
-template <int D, typename T>
-inline T& vector<D, T>::operator[](const unsigned int i)
-{
- return x[i];
-}
-
-
-// Dot
-// ~~~
-template <int D, typename T>
-inline double vector<D, T>::dot(const vector<D, T> &a) const
-{
- double d=0;
-
- for(int k=0; k<D; k++)
- d += x[k] * a.x[k];
-
- return d;
-}
-
-template <int D, typename T>
-inline double vector<D, T>::dot(const vector<D, T> &a, const vector<D, T> &b)
-{
- return a.dot(b);
-}
-
-
-// NormSquared
-// ~~~~~~~~~~~
-template <int D, typename T>
-inline double vector<D, T>::norm_squared() const
-{
- return dot(*this, *this);
-}
-
-template <int D, typename T>
-inline double vector<D, T>::norm_squared(const vector<D, T>& v)
-{
- return v.norm_squared();
-}
-
-
-// read
-// ~~~~
-template <int D, typename T>
-void vector<D, T>::read(std::ifstream& in)
-{
- in.read((char*)x, sizeof(T)*D);
-}
-
-// write
-// ~~~~~
-template <int D, typename T>
-void vector<D, T>::write(std::ofstream& out) const
-{
- out.write((const char*)x, sizeof(T)*D);
-}
-
-
-
-// Insertion
-// ~~~~~~~~~
-template <int D, typename T>
-std::ostream& operator<<(std::ostream& os, const vector<D, T>& v)
-{
- os << "(";
-
- for(int k=0; k<D-1; k++)
- os << v.x[k] << ", ";
-
- os << v.x[D-1] << ")";
-
- return os;
-}
-
-
-// ======================================================================
-// vector_field
-// ======================================================================
-
-// A field of V-vectors on a D dimensional manifold
-
-template<int V, int D, typename T=double>
-class vector_field {
-
- public:
- int elements;
-
- private:
- vector<V, T>* f;
- vector<D, int> size; // number of grid points for each dimension
- vector<D, int> offset;
-
- public:
-
- vector_field();
- vector_field(const vector<D, int>&);
- ~vector_field();
-
- vector<D, int> get_size() const;
- void set_size(const vector<D, int>&);
-
- vector<V, T>& get(const vector<D, int>&);
-
- void read(std::ifstream&);
- void write(std::ofstream&) const;
-
- static void swap(vector_field<V, D, T>&, vector_field<V, D, T>&);
-};
-
-
-// vector_field
-// ~~~~~~~~~~~~
-template<int V, int D, typename T>
-vector_field<V, D, T>::vector_field()
- : f(0), elements(0)
-{
-}
-
-
-// vector_field
-// ~~~~~~~~~~~~
-template<int V, int D, typename T>
-vector_field<V, D, T>::vector_field(const vector<D, int>& s)
- : f(0)
-{
- set_size(s);
-}
-
-// ~vector_field
-// ~~~~~~~~~~~~~
-template <int V, int D, typename T>
-vector_field<V, D, T>::~vector_field()
-{
- if(f != 0)
- delete[] f;
-}
-
-// get_size
-// ~~~~~~~~
-template<int V, int D, typename T>
-inline vector<D, int> vector_field<V, D, T>::get_size() const
-{
- return size;
-}
-
-// set_size
-// ~~~~~~~~
-template<int V, int D, typename T>
-void vector_field<V, D, T>::set_size(const vector<D, int>& s)
-{
- if(f != 0)
- delete[] f;
-
- size = s;
-
- elements = 1;
- for(int i=0; i<D; i++) {
- offset[i] = elements;
- elements *= size.x[i];
- }
-
- f = new vector<V, T>[elements];
-}
-
-// get
-// ~~~
-template<int V, int D, typename T>
-inline vector<V, T>& vector_field<V, D, T>::get(const vector<D, int>& pos)
-{
- int p=0;
- for(int i=0; i<D; i++)
- p += pos.x[i]*offset[i];
-
- return f[p];
-}
-
-
-// read
-// ~~~~
-template<int V, int D, typename T>
-void vector_field<V, D, T>::read(std::ifstream& in)
-{
- in.read((char*)f, elements*sizeof(T)*V);
-}
-
-
-// write
-// ~~~~~
-template<int V, int D, typename T>
-void vector_field<V, D, T>::write(std::ofstream& out) const
-{
- out.write((const char*)f, elements*sizeof(T)*V);
-}
-
-// swap
-// ~~~~
-template<int V, int D, typename T>
-void vector_field<V, D, T>::swap(vector_field<V, D, T>& v1,
- vector_field<V, D, T>& v2)
-{
- vector<V, T>* f;
-
- f = v1.f;
- v1.f = v2.f;
- v2.f = f;
-}
-
-
-
-#endif