From 901b9f16494f37890be17ef4bb66e6efc6873340 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Fri, 10 Feb 2017 12:18:11 -0500 Subject: changed code to rely on jst --- src/bound_set.c | 102 -------------------------------------------------------- 1 file changed, 102 deletions(-) delete mode 100644 src/bound_set.c (limited to 'src/bound_set.c') 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; -} -- cgit v1.2.3-70-g09d2