summaryrefslogtreecommitdiff
path: root/lib/src/graph.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'lib/src/graph.cpp')
-rw-r--r--lib/src/graph.cpp143
1 files changed, 75 insertions, 68 deletions
diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp
index a5063de..1475c40 100644
--- a/lib/src/graph.cpp
+++ b/lib/src/graph.cpp
@@ -34,23 +34,21 @@ graph::graph(unsigned Nx, unsigned Ny) {
vertices[i].r.x = (double)((1 + i / (Nx / 2)) % 2 + 2 * (i % (Nx / 2)));
vertices[i].r.y = (double)(i / (Nx / 2));
signed f = (1 + i / (Nx / 2)) % 2 == 1 ? 1 : -1;
- vertices[i].nd = {i, (i + Nx / 2) % nv, Nx / 2 * (i / (Nx / 2)) + ((i + f) % (Nx / 2)), (nv + i - Nx / 2) % nv};
- vertices[i].polygon = {
- {vertices[i].r.x - 1.0, vertices[i].r.y},
- {vertices[i].r.x, vertices[i].r.y - 1.0},
- {vertices[i].r.x + 1.0, vertices[i].r.y},
- {vertices[i].r.x, vertices[i].r.y + 1.0}
- };
+ vertices[i].nd = {i, (i + Nx / 2) % nv, Nx / 2 * (i / (Nx / 2)) + ((i + f) % (Nx / 2)),
+ (nv + i - Nx / 2) % nv};
+ vertices[i].polygon = {{vertices[i].r.x - 1.0, vertices[i].r.y},
+ {vertices[i].r.x, vertices[i].r.y - 1.0},
+ {vertices[i].r.x + 1.0, vertices[i].r.y},
+ {vertices[i].r.x, vertices[i].r.y + 1.0}};
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].nd = {i, (i + Nx / 2) % nv, Nx / 2 * (i / (Nx / 2)) + ((i + f) % (Nx / 2)), (nv + i - Nx / 2) % nv};
- dual_vertices[i].polygon = {
- {dual_vertices[i].r.x - 1.0, vertices[i].r.y},
- {dual_vertices[i].r.x, vertices[i].r.y - 1.0},
- {dual_vertices[i].r.x + 1.0, vertices[i].r.y},
- {dual_vertices[i].r.x, vertices[i].r.y + 1.0}
- };
+ dual_vertices[i].nd = {i, (i + Nx / 2) % nv, Nx / 2 * (i / (Nx / 2)) + ((i + f) % (Nx / 2)),
+ (nv + i - Nx / 2) % nv};
+ dual_vertices[i].polygon = {{dual_vertices[i].r.x - 1.0, vertices[i].r.y},
+ {dual_vertices[i].r.x, vertices[i].r.y - 1.0},
+ {dual_vertices[i].r.x + 1.0, vertices[i].r.y},
+ {dual_vertices[i].r.x, vertices[i].r.y + 1.0}};
}
for (unsigned y = 0; y < Ny; y++) {
@@ -66,7 +64,8 @@ graph::graph(unsigned Nx, unsigned Ny) {
unsigned dv1 = (Nx * y) / 2 + ((x + (y + 1) % 2) / 2) % (Nx / 2);
unsigned dv2 = ((Nx * (y + 1)) / 2 + ((x + y % 2) / 2) % (Nx / 2)) % nv;
- dual_edges.push_back({{dv1, dv2}, {0.5 + (double)x, 0.5 + (double)y}, {crossed_x, crossed_y}});
+ dual_edges.push_back(
+ {{dv1, dv2}, {0.5 + (double)x, 0.5 + (double)y}, {crossed_x, crossed_y}});
}
}
@@ -102,29 +101,23 @@ graph::graph(unsigned Nx, unsigned Ny) {
dual_vertices[vi].ne.push_back(i);
}
}
-
}
-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.";
+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;
-class clippingException: public std::exception
-{
- virtual const char* what() const throw()
- {
+class clippingException : public std::exception {
+ virtual const char* what() const throw() {
return "An interior site has a clipped edge and the tesselation is malformed.";
}
} clippingex;
-class triangleException: public std::exception
-{
- virtual const char* what() const throw()
- {
+class triangleException : public std::exception {
+ virtual const char* what() const throw() {
return "A dual-centered triangle has an impossible geometry and the tesselation is malformed.";
}
} triex;
@@ -147,29 +140,45 @@ unsigned get_triangle_signature(unsigned j1, unsigned j2, unsigned j3) {
if ((j1 == j2) && (j2 == j3)) {
return 0;
- } else if (((j1 == j2) && (x2 < x3) && (y2 == y3)) || ((j1 == j3) && (x3 < x2) && (y2 == y3)) || ((j2 == j3) && (x3 < x1) && (y1 == y3))) {
+ } else if (((j1 == j2) && (x2 < x3) && (y2 == y3)) || ((j1 == j3) && (x3 < x2) && (y2 == y3)) ||
+ ((j2 == j3) && (x3 < x1) && (y1 == y3))) {
return 1;
- } else if (((j1 == j2) && (x2 > x3) && (y2 == y3)) || ((j1 == j3) && (x3 > x2) && (y2 == y3)) || ((j2 == j3) && (x3 > x1) && (y1 == y3))) {
+ } else if (((j1 == j2) && (x2 > x3) && (y2 == y3)) || ((j1 == j3) && (x3 > x2) && (y2 == y3)) ||
+ ((j2 == j3) && (x3 > x1) && (y1 == y3))) {
return 2;
- } else if (((j1 == j2) && (y2 < y3) && (x2 == x3)) || ((j1 == j3) && (y3 < y2) && (x2 == x3)) || ((j2 == j3) && (y3 < y1) && (x1 == x3))) {
+ } else if (((j1 == j2) && (y2 < y3) && (x2 == x3)) || ((j1 == j3) && (y3 < y2) && (x2 == x3)) ||
+ ((j2 == j3) && (y3 < y1) && (x1 == x3))) {
return 3;
- } else if (((j1 == j2) && (y2 > y3) && (x2 == x3)) || ((j1 == j3) && (y3 > y2) && (x2 == x3)) || ((j2 == j3) && (y3 > y1) && (x1 == x3))) {
+ } else if (((j1 == j2) && (y2 > y3) && (x2 == x3)) || ((j1 == j3) && (y3 > y2) && (x2 == x3)) ||
+ ((j2 == j3) && (y3 > y1) && (x1 == x3))) {
return 4;
- } else if (((j1 == j2) && (x2 < x3) && (y2 < y3)) || ((j1 == j3) && (x3 < x2) && (y3 < y2)) || ((j2 == j3) && (x3 < x1) && (y3 < y1))) {
+ } else if (((j1 == j2) && (x2 < x3) && (y2 < y3)) || ((j1 == j3) && (x3 < x2) && (y3 < y2)) ||
+ ((j2 == j3) && (x3 < x1) && (y3 < y1))) {
return 5;
- } else if (((j1 == j2) && (x2 < x3) && (y2 > y3)) || ((j1 == j3) && (x3 < x2) && (y3 > y2)) || ((j2 == j3) && (x3 < x1) && (y3 > y1))) {
+ } else if (((j1 == j2) && (x2 < x3) && (y2 > y3)) || ((j1 == j3) && (x3 < x2) && (y3 > y2)) ||
+ ((j2 == j3) && (x3 < x1) && (y3 > y1))) {
return 6;
- } else if (((j1 == j2) && (x2 > x3) && (y2 < y3)) || ((j1 == j3) && (x3 > x2) && (y3 < y2)) || ((j2 == j3) && (x3 > x1) && (y3 < y1))) {
+ } else if (((j1 == j2) && (x2 > x3) && (y2 < y3)) || ((j1 == j3) && (x3 > x2) && (y3 < y2)) ||
+ ((j2 == j3) && (x3 > x1) && (y3 < y1))) {
return 7;
- } else if (((j1 == j2) && (x2 > x3) && (y2 > y3)) || ((j1 == j3) && (x3 > x2) && (y3 > y2)) || ((j2 == j3) && (x3 > x1) && (y3 > y1))) {
+ } else if (((j1 == j2) && (x2 > x3) && (y2 > y3)) || ((j1 == j3) && (x3 > x2) && (y3 > y2)) ||
+ ((j2 == j3) && (x3 > x1) && (y3 > y1))) {
return 8;
- } else if (((x1 == x2) && (x2 < x3) && ((y1 < y3) || (y2 < y3))) || ((x1 == x3) && (x3 < x2) && ((y1 < y2) || (y3 < y2))) || ((x2 == x3) && (x2 < x1) && ((y2 < y1) || (y3 < y1)))) {
+ } else if (((x1 == x2) && (x2 < x3) && ((y1 < y3) || (y2 < y3))) ||
+ ((x1 == x3) && (x3 < x2) && ((y1 < y2) || (y3 < y2))) ||
+ ((x2 == x3) && (x2 < x1) && ((y2 < y1) || (y3 < y1)))) {
return 9;
- } else if (((x1 == x2) && (x2 > x3) && ((y1 < y3) || (y2 < y3))) || ((x1 == x3) && (x3 > x2) && ((y1 < y2) || (y3 < y2))) || ((x2 == x3) && (x2 > x1) && ((y2 < y1) || (y3 < y1)))) {
+ } else if (((x1 == x2) && (x2 > x3) && ((y1 < y3) || (y2 < y3))) ||
+ ((x1 == x3) && (x3 > x2) && ((y1 < y2) || (y3 < y2))) ||
+ ((x2 == x3) && (x2 > x1) && ((y2 < y1) || (y3 < y1)))) {
return 10;
- } else if (((x1 == x2) && (x2 < x3) && ((y1 > y3) || (y2 > y3))) || ((x1 == x3) && (x3 < x2) && ((y1 > y2) || (y3 > y2))) || ((x2 == x3) && (x2 < x1) && ((y2 > y1) || (y3 > y1)))) {
+ } else if (((x1 == x2) && (x2 < x3) && ((y1 > y3) || (y2 > y3))) ||
+ ((x1 == x3) && (x3 < x2) && ((y1 > y2) || (y3 > y2))) ||
+ ((x2 == x3) && (x2 < x1) && ((y2 > y1) || (y3 > y1)))) {
return 11;
- } else if (((x1 == x2) && (x2 > x3) && ((y1 > y3) || (y2 > y3))) || ((x1 == x3) && (x3 > x2) && ((y1 > y2) || (y3 > y2))) || ((x2 == x3) && (x2 > x1) && ((y2 > y1) || (y3 > y1)))) {
+ } else if (((x1 == x2) && (x2 > x3) && ((y1 > y3) || (y2 > y3))) ||
+ ((x1 == x3) && (x3 > x2) && ((y1 > y2) || (y3 > y2))) ||
+ ((x2 == x3) && (x2 > x1) && ((y2 > y1) || (y3 > y1)))) {
return 12;
} else {
throw triex;
@@ -200,7 +209,7 @@ void graph::helper(unsigned nv, std::mt19937& rng) {
// the coordinates of the lattice, from which a delaunay triangulation
// and corresponding voronoi tessellation will be built. Everyone is in the
// rectangle (0, 0) (Lx, Ly)
- for (vertex &v : vertices) {
+ for (vertex& v : vertices) {
v = {{L.x * d(rng), L.y * d(rng)}};
}
@@ -243,7 +252,7 @@ void graph::helper(unsigned nv, std::mt19937& rng) {
ep = ep->next;
}
// for each edge on the site's cell
- while(e) {
+ while (e) {
// 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
@@ -275,7 +284,8 @@ void graph::helper(unsigned nv, std::mt19937& rng) {
if (it1 == known_vertices.end()) {
vi1 = dual_vertices.size();
- dual_vertices.push_back({{mod(e->pos[0].x, L.x), mod(e->pos[0].y, L.y)},{i1},{},{{site->p.x, site->p.y}}});
+ dual_vertices.push_back(
+ {{mod(e->pos[0].x, L.x), mod(e->pos[0].y, L.y)}, {i1}, {}, {{site->p.x, site->p.y}}});
known_vertices[t1] = vi1;
} else {
vi1 = it1->second;
@@ -284,7 +294,7 @@ void graph::helper(unsigned nv, std::mt19937& rng) {
}
vertices[i1].nd.push_back(vi1);
-
+
if (i1 < i2 || (i1 == i2 && !self_bonded)) {
if (i1 == i2) {
self_bonded = true;
@@ -292,12 +302,11 @@ void graph::helper(unsigned nv, std::mt19937& rng) {
bool crossed_x = x2 != 1;
bool crossed_y = y2 != 1;
edges.push_back({{i1, i2},
- {mod((site->p.x + neighbor->p.x) / 2, L.x),
- mod((site->p.y + neighbor->p.y) / 2, L.y)},
- {crossed_x, crossed_y}
- });
+ {mod((site->p.x + neighbor->p.x) / 2, L.x),
+ mod((site->p.y + neighbor->p.y) / 2, L.y)},
+ {crossed_x, crossed_y}});
- jcv_graphedge *en;
+ jcv_graphedge* en;
if (e->next == NULL) {
en = site->edges;
} else {
@@ -322,20 +331,21 @@ void graph::helper(unsigned nv, std::mt19937& rng) {
if (it2 == known_vertices.end()) {
vi2 = dual_vertices.size();
- dual_vertices.push_back({{mod(e->pos[1].x, L.x), mod(e->pos[1].y, L.y)},{}});
+ dual_vertices.push_back({{mod(e->pos[1].x, L.x), mod(e->pos[1].y, L.y)}, {}});
known_vertices[t2] = vi2;
} else {
vi2 = it2->second;
}
- bool dcrossed_x = (unsigned)floor(e->pos[0].x / L.x) != (unsigned)floor(e->pos[1].x / L.x);
- bool dcrossed_y = (unsigned)floor(e->pos[0].y / L.y) != (unsigned)floor(e->pos[1].y / L.y);
+ bool dcrossed_x =
+ (unsigned)floor(e->pos[0].x / L.x) != (unsigned)floor(e->pos[1].x / L.x);
+ bool dcrossed_y =
+ (unsigned)floor(e->pos[0].y / L.y) != (unsigned)floor(e->pos[1].y / L.y);
dual_edges.push_back({{vi1, vi2},
- {mod((e->pos[0].x + e->pos[1].x) / 2, L.x),
- mod((e->pos[0].y + e->pos[1].y) / 2, L.y)},
- {dcrossed_x, dcrossed_y}
- });
+ {mod((e->pos[0].x + e->pos[1].x) / 2, L.x),
+ mod((e->pos[0].y + e->pos[1].y) / 2, L.y)},
+ {dcrossed_x, dcrossed_y}});
}
ep = e;
@@ -348,7 +358,7 @@ void graph::helper(unsigned nv, std::mt19937& rng) {
throw eulerex;
}
- for (vertex &v : dual_vertices) {
+ 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;
@@ -434,33 +444,31 @@ void graph::helper(unsigned nv, std::mt19937& rng) {
jcv_diagram_free(&diagram);
}
-graph::coordinate reverse(const graph::coordinate &x) {
- return {x.y, x.x};
-}
+graph::coordinate reverse(const graph::coordinate& x) { return {x.y, x.x}; }
graph const graph::rotate() {
graph g2(*this);
- for (graph::vertex &v : g2.vertices) {
+ for (graph::vertex& v : g2.vertices) {
v.r = reverse(v.r);
- for (graph::coordinate &r : v.polygon) {
+ for (graph::coordinate& r : v.polygon) {
r = reverse(r);
}
}
- for (graph::edge &e : g2.edges) {
+ for (graph::edge& e : g2.edges) {
e.r = reverse(e.r);
e.crossings = {e.crossings[1], e.crossings[0]};
}
- for (graph::vertex &v : g2.dual_vertices) {
+ for (graph::vertex& v : g2.dual_vertices) {
v.r = reverse(v.r);
- for (graph::coordinate &r : v.polygon) {
+ for (graph::coordinate& r : v.polygon) {
r = reverse(r);
}
}
- for (graph::edge &e : g2.dual_edges) {
+ for (graph::edge& e : g2.dual_edges) {
e.r = reverse(e.r);
e.crossings = {e.crossings[1], e.crossings[0]};
}
@@ -469,4 +477,3 @@ graph const graph::rotate() {
return g2;
}
-