diff options
Diffstat (limited to 'lib/bound_set.c')
| -rw-r--r-- | lib/bound_set.c | 178 |
1 files changed, 93 insertions, 85 deletions
diff --git a/lib/bound_set.c b/lib/bound_set.c index 65f1e6f..32d0fb7 100644 --- a/lib/bound_set.c +++ b/lib/bound_set.c @@ -2,101 +2,109 @@ #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; - + 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))); + 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); + 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); - } + 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; + 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; +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; } |