From 09200a607661f739782a966807d31345485e2c41 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Thu, 20 Dec 2018 12:20:19 -0500 Subject: added animation example, and did many fixes to the voronoi system --- lib/src/graph.cpp | 217 +++++++++++++++++++++++++++++++++++++++------------- lib/src/network.cpp | 28 +++---- 2 files changed, 178 insertions(+), 67 deletions(-) (limited to 'lib/src') diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp index f25cf87..10f62b7 100644 --- a/lib/src/graph.cpp +++ b/lib/src/graph.cpp @@ -1,7 +1,6 @@ #include "graph.hpp" #include -#include #define JC_VORONOI_IMPLEMENTATION #define JCV_REAL_TYPE double @@ -24,9 +23,21 @@ graph::graph(unsigned int Nx, unsigned int Ny) { for (unsigned int i = 0; i < nv; i++) { vertices[i].r.x = (double)((1 + i / (Nx / 2)) % 2 + 2 * (i % (Nx / 2))); vertices[i].r.y = (double)(i / (Nx / 2)); + vertices[i].polygon = { + {vertices[i].r.x - 0.5, vertices[i].r.y}, + {vertices[i].r.x, vertices[i].r.y - 0.5}, + {vertices[i].r.x + 0.5, vertices[i].r.y}, + {vertices[i].r.x, vertices[i].r.y + 0.5} + }; dual_vertices[i].r.x = (double)((i / (Nx / 2)) % 2 + 2 * (i % (Nx / 2))); dual_vertices[i].r.y = (double)(i / (Nx / 2)); + dual_vertices[i].polygon = { + {dual_vertices[i].r.x - 0.5, vertices[i].r.y}, + {dual_vertices[i].r.x, vertices[i].r.y - 0.5}, + {dual_vertices[i].r.x + 0.5, vertices[i].r.y}, + {dual_vertices[i].r.x, vertices[i].r.y + 0.5} + }; } for (unsigned int x = 0; x < Ny; x++) { @@ -34,40 +45,70 @@ graph::graph(unsigned int Nx, unsigned int Ny) { unsigned int v1 = (Nx * x) / 2 + ((y + x % 2) / 2) % (Nx / 2); unsigned int v2 = ((Nx * (x + 1)) / 2 + ((y + (x + 1) % 2) / 2) % (Nx / 2)) % nv; - edges.push_back({v1, v2}); + edges.push_back({{v1, v2}, {0.5 + (double)y, 0.5 + (double)x}}); unsigned int dv1 = (Nx * x) / 2 + ((y + (x + 1) % 2) / 2) % (Nx / 2); unsigned int dv2 = ((Nx * (x + 1)) / 2 + ((y + x % 2) / 2) % (Nx / 2)) % nv; - dual_edges.push_back({dv1, dv2}); + dual_edges.push_back({{dv1, dv2}, {0.5 + (double)y, 0.5 + (double)x}}); } } } +namespace std +{ + template + struct hash > + { + typedef array argument_type; + typedef size_t result_type; + + result_type operator()(const argument_type& a) const + { + hash hasher; + result_type h = 0; + for (result_type i = 0; i < N; ++i) + { + h = h * 31 + hasher(a[i]); + } + return h; + } + }; +} + +class eulerException: public std::exception +{ + virtual const char* what() const throw() + { + return "The voronoi tessellation generated has the incorrect Euler characteristic for a torus and is malformed."; + } +} eulerex; + + graph::graph(unsigned int Nx, unsigned int Ny, std::mt19937& rng, double spread) { L = {(double)Nx, (double)Ny}; - unsigned int dnv = Nx * Ny / 2; + unsigned int nv = Nx * Ny / 2; - dual_vertices.resize(dnv); + vertices.resize(nv); std::normal_distribution d(0.0, spread); - // the coordinates of the dual lattice, from which a delaunay triangulation - // and corresponding voronoi tesselation will be built - for (unsigned int i = 0; i < dnv; i++) { + // the coordinates of the lattice, from which a delaunay triangulation + // and corresponding voronoi tessellation will be built + for (unsigned int i = 0; i < nv; i++) { double rx = (double)((i / (Nx / 2)) % 2 + 2 * (i % (Nx / 2))) + d(rng); double ry = (double)(i / (Nx / 2)) + d(rng); - dual_vertices[i] = {{fmod(L.x + rx, L.x), fmod(L.y + ry, L.y)}}; + vertices[i] = {{fmod(L.x + rx, L.x), fmod(L.y + ry, L.y)}}; } - // to make the resulting tesselation periodic, we tile eight (!) copies of - // the original dual points for a total of nine. note that the index of each + // to make the resulting tessellation periodic, we tile eight (!) copies of + // the original points for a total of nine. note that the index of each // point quotient 9 is equal to the original index (we'll use this to // translate back) std::vector points; - points.reserve(9 * dnv); - for (const vertex &v : dual_vertices) { + points.reserve(9 * nv); + for (const vertex &v : vertices) { points.push_back({v.r.x, v.r.y}); points.push_back({v.r.x + L.x, v.r.y + 0.0}); points.push_back({v.r.x - L.x, v.r.y + 0.0}); @@ -82,81 +123,151 @@ graph::graph(unsigned int Nx, unsigned int Ny, std::mt19937& rng, double spread) jcv_diagram diagram; memset(&diagram, 0, sizeof(jcv_diagram)); - jcv_diagram_generate(9 * dnv, points.data(), NULL, &diagram); + jcv_diagram_generate(9 * nv, points.data(), NULL, &diagram); const jcv_site* sites = jcv_diagram_get_sites(&diagram); - struct coorcomp { - bool operator() (const coordinate& lhs, const coordinate& rhs) const - {return lhs.x& lhs, const std::array& rhs) const + {if (lhs[0] != lhs[1]) return lhs[0] < lhs[1]; + else return rhs[0] < rhs[1]; + } }; - std::map known_vertices; + std::unordered_map, unsigned int> known_vertices; for (int i = 0; i < diagram.numsites; i++) { const jcv_site* site = &sites[i]; // we only care about processing the cells of our original, central sites if (site->index % 9 == 0) { + unsigned int i1 = (unsigned int)(site->index / 9); const jcv_graphedge* e = site->edges; + const jcv_graphedge* ep = site->edges; + while (ep->next) { + ep = ep->next; + } // for each edge on the site's cell while(e) { - // assess whether the dual edge needs to be added. only neighboring - // dual sites whose index is larger than ours are given edges, so we + // assess whether the edge needs to be added. only neighboring + // sites whose index is larger than ours are given edges, so we // don't double up later const jcv_site* neighbor = e->neighbor; - unsigned int index = (unsigned int)(site->index / 9); - unsigned int real_index = (unsigned int)(neighbor->index / 9); - - if (index < real_index) { - dual_edges.push_back({index, real_index}); - coordinate r1 = {1e-10 * round(1e10 * fmod(L.x + e->pos[0].x, L.x)), 1e-10 * round(1e10 * fmod(L.y + e->pos[0].y, L.y))}; - coordinate r2 = {1e-10 * round(1e10 * fmod(L.x + e->pos[1].x, L.x)), 1e-10 * round(1e10 * fmod(L.y + e->pos[1].y, L.y))}; + unsigned int i2 = (unsigned int)(neighbor->index / 9); + + vertices[i1].polygon.push_back({e->pos[0].x, e->pos[0].y}); + + unsigned int i3p = (unsigned int)(ep->neighbor->index / 9); + std::array t1 = {i1, i2, i3p}; + std::sort(t1.begin(), t1.end()); - std::map::iterator it1 = known_vertices.find(r1); - std::map::iterator it2 = known_vertices.find(r2); + auto it1 = known_vertices.find(t1); - unsigned int vi1, vi2; + unsigned int vi1; - if (it1 == known_vertices.end()) { - vi1 = vertices.size(); - vertices.push_back({r1}); - known_vertices[r1] = vi1; + if (it1 == known_vertices.end()) { + vi1 = dual_vertices.size(); + dual_vertices.push_back({{fmod(L.x + e->pos[0].x, L.x), fmod(L.y + e->pos[0].y, L.y)},{{site->p.x, site->p.y}}}); + known_vertices[t1] = vi1; + } else { + vi1 = it1->second; + dual_vertices[vi1].polygon.push_back({site->p.x, site->p.y}); + } + + if (i1 < i2) { + edges.push_back({{i1, i2}, + {fmod(L.x + (site->p.x + neighbor->p.x) / 2, L.x), + fmod(L.y + (site->p.y + neighbor->p.y) / 2, L.y)} + }); + + jcv_graphedge *en; + if (e->next == NULL) { + en = site->edges; } else { - vi1 = it1->second; + en = e->next; } + + unsigned int i3n = (unsigned int)(en->neighbor->index / 9); + std::array t2 = {i1, i2, i3n}; + std::sort(t2.begin(), t2.end()); + + auto it2 = known_vertices.find(t2); + + unsigned int vi2; + if (it2 == known_vertices.end()) { - vi2 = vertices.size(); - vertices.push_back({r2}); - known_vertices[r2] = vi2; + vi2 = dual_vertices.size(); + dual_vertices.push_back({{fmod(L.x + e->pos[1].x, L.x), fmod(L.y + e->pos[1].y, L.y)},{}}); + known_vertices[t2] = vi2; } else { vi2 = it2->second; } - edges.push_back({vi1, vi2}); + dual_edges.push_back({{vi1, vi2}, + {fmod(L.x + (e->pos[0].x + e->pos[1].x) / 2, L.x), + fmod(L.y + (e->pos[0].y + e->pos[1].y) / 2, L.y)} + }); } + ep = e; e = e->next; } } } - for (edge &e : edges) { - std::cout << e[0] << " " << e[1] << " "; - } - std::cout << "\n"; - for (vertex &v : vertices) { - std::cout << v.r.x << " " << v.r.y << " "; + if (edges.size() != vertices.size() + dual_vertices.size()) { + throw eulerex; } - std::cout << "\n"; - for (edge &e : dual_edges) { - std::cout << e[0] << " " << e[1] << " "; - } - std::cout << "\n"; - for (vertex &v : dual_vertices) { - std::cout << v.r.x << " " << v.r.y << " "; + + for (vertex &v : dual_vertices) { + if (fabs(v.polygon[0].x - v.polygon[1].x) > L.x / 2) { + if (v.polygon[0].x < L.x / 2) { + v.polygon[0].x += L.x; + } else { + v.polygon[1].x += L.x; + } + } + + if (fabs(v.polygon[0].y - v.polygon[1].y) > L.y / 2) { + if (v.polygon[0].y < L.y / 2) { + v.polygon[0].y += L.y; + } else { + v.polygon[1].y += L.y; + } + } + + if (fabs(v.polygon[2].x - v.polygon[1].x) > L.x / 2) { + if (v.polygon[2].x < L.x / 2) { + v.polygon[2].x += L.x; + } else { + v.polygon[1].x += L.x; + } + } + + if (fabs(v.polygon[2].y - v.polygon[1].y) > L.y / 2) { + if (v.polygon[2].y < L.y / 2) { + v.polygon[2].y += L.y; + } else { + v.polygon[1].y += L.y; + } + } + + if (fabs(v.polygon[2].x - v.polygon[0].x) > L.x / 2) { + if (v.polygon[2].x < L.x / 2) { + v.polygon[2].x += L.x; + } else { + v.polygon[0].x += L.x; + } + } + + if (fabs(v.polygon[2].y - v.polygon[0].y) > L.y / 2) { + if (v.polygon[2].y < L.y / 2) { + v.polygon[2].y += L.y; + } else { + v.polygon[0].y += L.y; + } + } } - std::cout << "\n"; jcv_diagram_free(&diagram); } diff --git a/lib/src/network.cpp b/lib/src/network.cpp index 9cb1007..4d4ed2d 100644 --- a/lib/src/network.cpp +++ b/lib/src/network.cpp @@ -4,14 +4,14 @@ network::network(const graph& G, cholmod_common *c) : c(c), G(G), fuses(G.edges.size(), false), thresholds(G.edges.size(), 1) { b = CHOL_F(zeros)(G.vertices.size(), 1, CHOLMOD_REAL, c); for (unsigned int i = 0; i < G.edges.size(); i++) { - double v0y = G.vertices[G.edges[i][0]].r.y; - double v1y = G.vertices[G.edges[i][1]].r.y; + double v0y = G.vertices[G.edges[i].v[0]].r.y; + double v1y = G.vertices[G.edges[i].v[1]].r.y; if (fabs(v0y - v1y) > G.L.y / 2) { bool ind = v0y < v1y ? 0 : 1; - ((double *)b->x)[G.edges[i][ind]] += 1.0; - ((double *)b->x)[G.edges[i][!ind]] -= 1.0; + ((double *)b->x)[G.edges[i].v[ind]] += 1.0; + ((double *)b->x)[G.edges[i].v[!ind]] -= 1.0; } } @@ -30,8 +30,8 @@ network::network(const graph& G, cholmod_common *c) : c(c), G(G), fuses(G.edges. unsigned int terms = G.vertices.size(); for (unsigned int i = 0; i < G.edges.size(); i++) { - unsigned int v0 = G.edges[i][0]; - unsigned int v1 = G.edges[i][1]; + unsigned int v0 = G.edges[i].v[0]; + unsigned int v1 = G.edges[i].v[1]; unsigned int s0 = v0 < v1 ? v0 : v1; unsigned int s1 = v0 < v1 ? v1 : v0; @@ -60,11 +60,11 @@ network::network(const graph& G, cholmod_common *c) : c(c), G(G), fuses(G.edges. for (unsigned int i = 0; i < G.edges.size(); i++) { ((CHOL_INT *)t->i)[2 * i] = i; - ((CHOL_INT *)t->j)[2 * i] = G.edges[i][0]; + ((CHOL_INT *)t->j)[2 * i] = G.edges[i].v[0]; ((double *)t->x)[2 * i] = 1.0; ((CHOL_INT *)t->i)[2 * i + 1] = i; - ((CHOL_INT *)t->j)[2 * i + 1] = G.edges[i][1]; + ((CHOL_INT *)t->j)[2 * i + 1] = G.edges[i].v[1]; ((double *)t->x)[2 * i + 1] = -1.0; } @@ -99,8 +99,8 @@ void network::set_thresholds(double beta, std::mt19937& rng) { void network::break_edge(unsigned int e, bool unbreak) { fuses[e] = !unbreak; - unsigned int v0 = G.edges[e][0]; - unsigned int v1 = G.edges[e][1]; + unsigned int v0 = G.edges[e].v[0]; + unsigned int v1 = G.edges[e].v[1]; unsigned int n = factor->n; @@ -141,8 +141,8 @@ void network::break_edge(unsigned int e, bool unbreak) { if (fabs(v0y - v1y) > G.L.y / 2) { bool ind = v0y < v1y ? unbreak : !unbreak; - ((double *)b->x)[G.edges[e][ind]] -= 1.0; - ((double *)b->x)[G.edges[e][!ind]] += 1.0; + ((double *)b->x)[G.edges[e].v[ind]] -= 1.0; + ((double *)b->x)[G.edges[e].v[!ind]] += 1.0; } } @@ -169,8 +169,8 @@ current_info network::get_current_info() { } else { currents[i] = ((double *)y->x)[i]; - double v0y = G.vertices[G.edges[i][0]].r.y; - double v1y = G.vertices[G.edges[i][1]].r.y; + double v0y = G.vertices[G.edges[i].v[0]].r.y; + double v1y = G.vertices[G.edges[i].v[1]].r.y; if (fabs(v0y - v1y) > G.L.y / 2) { if (v0y > v1y) { -- cgit v1.2.3-70-g09d2