diff options
Diffstat (limited to 'src')
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 |