summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--lib/include/array_hash.hpp27
-rw-r--r--lib/include/graph.hpp2
-rw-r--r--lib/include/network.hpp1
-rw-r--r--lib/src/graph.cpp34
-rw-r--r--lib/src/network.cpp33
-rw-r--r--src/analysis_tools.cpp6
-rw-r--r--src/animate.cpp48
7 files changed, 87 insertions, 64 deletions
diff --git a/lib/include/array_hash.hpp b/lib/include/array_hash.hpp
new file mode 100644
index 0000000..3b35e4a
--- /dev/null
+++ b/lib/include/array_hash.hpp
@@ -0,0 +1,27 @@
+
+#pragma once
+
+#include <array>
+#include <unordered_map>
+
+namespace std
+{
+ template<typename T, size_t N>
+ struct hash<array<T, N> >
+ {
+ typedef array<T, N> argument_type;
+ typedef size_t result_type;
+
+ result_type operator()(const argument_type& a) const
+ {
+ hash<T> hasher;
+ result_type h = 0;
+ for (result_type i = 0; i < N; ++i)
+ {
+ h = h * 31 + hasher(a[i]);
+ }
+ return h;
+ }
+ };
+}
+
diff --git a/lib/include/graph.hpp b/lib/include/graph.hpp
index c1ad4f0..ea08214 100644
--- a/lib/include/graph.hpp
+++ b/lib/include/graph.hpp
@@ -9,6 +9,7 @@
#include <list>
#include <exception>
+#include "array_hash.hpp"
class graph {
public:
@@ -25,6 +26,7 @@ class graph {
typedef struct edge {
std::array<unsigned int, 2> v;
coordinate r;
+ std::array<bool, 2> crossings;
} edge;
coordinate L;
diff --git a/lib/include/network.hpp b/lib/include/network.hpp
index d955e43..70173e6 100644
--- a/lib/include/network.hpp
+++ b/lib/include/network.hpp
@@ -12,6 +12,7 @@
#include "graph.hpp"
#include "hooks.hpp"
#include "current_info.hpp"
+#include "array_hash.hpp"
#ifdef FRACTURE_LONGINT
diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp
index 2fee5fd..7564fa7 100644
--- a/lib/src/graph.cpp
+++ b/lib/src/graph.cpp
@@ -6,7 +6,6 @@
#define JCV_REAL_TYPE double
#define JCV_ATAN2 atan2
#define JCV_FLT_MAX 1.7976931348623157E+308
-#define NDEBUG
#include <jc_voronoi.h>
// actual mod for floats
@@ -65,27 +64,6 @@ graph::graph(unsigned int Nx, unsigned int Ny) {
}
}
-namespace std
-{
- template<typename T, size_t N>
- struct hash<array<T, N> >
- {
- typedef array<T, N> argument_type;
- typedef size_t result_type;
-
- result_type operator()(const argument_type& a) const
- {
- hash<T> 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()
@@ -248,9 +226,13 @@ graph::graph(double Lx, double Ly, std::mt19937& rng, double relax, double step)
}
if (i1 < i2) {
+ unsigned int nn = neighbor->index % 9;
+ bool crossed_x = (nn == 1) || (nn == 2) || (nn > 4);
+ bool crossed_y = nn > 2;
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)}
+ mod((site->p.y + neighbor->p.y) / 2, L.y)},
+ {crossed_x, crossed_y}
});
jcv_graphedge *en;
@@ -276,9 +258,13 @@ graph::graph(double Lx, double Ly, std::mt19937& rng, double relax, double step)
vi2 = it2->second;
}
+ bool dcrossed_x = (unsigned int)floor(e->pos[0].x / L.x) != (unsigned int)floor(e->pos[1].x / L.x);
+ bool dcrossed_y = (unsigned int)floor(e->pos[0].y / L.y) != (unsigned int)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)}
+ mod((e->pos[0].y + e->pos[1].y) / 2, L.y)},
+ {dcrossed_x, dcrossed_y}
});
}
diff --git a/lib/src/network.cpp b/lib/src/network.cpp
index 4d4ed2d..3812f43 100644
--- a/lib/src/network.cpp
+++ b/lib/src/network.cpp
@@ -7,7 +7,7 @@ network::network(const graph& G, cholmod_common *c) : c(c), G(G), fuses(G.edges.
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 (G.edges[i].crossings[1]) {
bool ind = v0y < v1y ? 0 : 1;
((double *)b->x)[G.edges[i].v[ind]] += 1.0;
@@ -19,8 +19,6 @@ network::network(const graph& G, cholmod_common *c) : c(c), G(G), fuses(G.edges.
cholmod_triplet *t = CHOL_F(allocate_triplet)(G.vertices.size(), G.vertices.size(), nnz, 1, CHOLMOD_REAL, c);
- t->nnz = nnz;
-
for (unsigned int i = 0; i < G.vertices.size(); i++) {
((CHOL_INT *)t->i)[i] = i;
((CHOL_INT *)t->j)[i] = i;
@@ -29,26 +27,37 @@ network::network(const graph& G, cholmod_common *c) : c(c), G(G), fuses(G.edges.
unsigned int terms = G.vertices.size();
+ std::unordered_map<std::array<unsigned int, 2>, unsigned int> known_edges;
+
for (unsigned int i = 0; i < G.edges.size(); i++) {
unsigned int v0 = G.edges[i].v[0];
unsigned int v1 = G.edges[i].v[1];
+ ((double *)t->x)[v0]++;
+ ((double *)t->x)[v1]++;
+
unsigned int s0 = v0 < v1 ? v0 : v1;
unsigned int s1 = v0 < v1 ? v1 : v0;
- ((CHOL_INT *)t->i)[terms] = s1;
- ((CHOL_INT *)t->j)[terms] = s0;
- ((double *)t->x)[terms] = -1.0;
+ auto it = known_edges.find({s0, s1});
- ((double *)t->x)[v0]++;
- ((double *)t->x)[v1]++;
+ if (it == known_edges.end()) {
+ ((CHOL_INT *)t->i)[terms] = s1;
+ ((CHOL_INT *)t->j)[terms] = s0;
+ ((double *)t->x)[terms] = -1.0;
- terms++;
+ known_edges[{s0, s1}] = terms;
+ terms++;
+ } else {
+ ((double *)t->x)[it->second] -= 1.0;
+ }
}
((double *)t->x)[0]++;
- cholmod_sparse *laplacian = CHOL_F(triplet_to_sparse)(t, nnz, c);
+ t->nnz = terms;
+
+ cholmod_sparse *laplacian = CHOL_F(triplet_to_sparse)(t, terms, c);
CHOL_F(free_triplet)(&t, c);
factor = CHOL_F(analyze)(laplacian, c);
CHOL_F(factorize)(laplacian, factor, c);
@@ -138,7 +147,7 @@ void network::break_edge(unsigned int e, bool unbreak) {
double v0y = G.vertices[v0].r.y;
double v1y = G.vertices[v1].r.y;
- if (fabs(v0y - v1y) > G.L.y / 2) {
+ if (G.edges[e].crossings[1]) {
bool ind = v0y < v1y ? unbreak : !unbreak;
((double *)b->x)[G.edges[e].v[ind]] -= 1.0;
@@ -172,7 +181,7 @@ current_info network::get_current_info() {
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 (G.edges[i].crossings[1]) {
if (v0y > v1y) {
currents[i] += 1.0;
total_current += currents[i];
diff --git a/src/analysis_tools.cpp b/src/analysis_tools.cpp
index 1e799bb..34e4ea0 100644
--- a/src/analysis_tools.cpp
+++ b/src/analysis_tools.cpp
@@ -102,12 +102,10 @@ std::list<unsigned int> find_minimal_crack(const Graph& G, const network& n) {
std::array<unsigned int, 2> crossing_count{0,0};
for (auto edge : cycle) {
- double dx = fabs(n.G.dual_vertices[n.G.dual_edges[edge].v[0]].r.x - n.G.dual_vertices[n.G.dual_edges[edge].v[1]].r.x);
- if (dx > n.G.L.x / 2) {
+ if (n.G.dual_edges[edge].crossings[0]) {
crossing_count[0]++;
}
- double dy = fabs(n.G.dual_vertices[n.G.dual_edges[edge].v[0]].r.y - n.G.dual_vertices[n.G.dual_edges[edge].v[1]].r.y);
- if (dy > n.G.L.y / 2) {
+ if (n.G.dual_edges[edge].crossings[1]) {
crossing_count[1]++;
}
}
diff --git a/src/animate.cpp b/src/animate.cpp
index 26f4996..a24f0f2 100644
--- a/src/animate.cpp
+++ b/src/animate.cpp
@@ -25,15 +25,15 @@ void animate::pre_fracture(const network &n) {
graph::coordinate r1 = n.G.vertices[n.G.edges[i].v[0]].r;
graph::coordinate r2 = n.G.vertices[n.G.edges[i].v[1]].r;
- if (fabs(r1.x - r2.x) > n.G.L.x / 2) {
- if (r1.x < n.G.L.x / 2) {
+ if (n.G.edges[i].crossings[0]) {
+ if (r1.x < r2.x) {
r1.x += n.G.L.x;
} else {
r2.x += n.G.L.x;
}
}
- if (fabs(r1.y - r2.y) > n.G.L.y / 2) {
- if (r1.y < n.G.L.y / 2) {
+ if (n.G.edges[i].crossings[1]) {
+ if (r1.y < r2.y) {
r1.y += n.G.L.y;
} else {
r2.y += n.G.L.y;
@@ -63,16 +63,16 @@ void animate::bond_broken(const network& n, const current_info& cur, unsigned in
graph::coordinate r1 = n.G.vertices[n.G.edges[i].v[0]].r;
graph::coordinate r2 = n.G.vertices[n.G.edges[i].v[1]].r;
- if (fabs(r1.x - r2.x) > n.G.L.x / 2) {
- if (r1.x < n.G.L.x / 2) {
- r2.x -= n.G.L.x;
+ if (n.G.edges[i].crossings[0]) {
+ if (r1.x < r2.x) {
+ r1.x += n.G.L.x;
} else {
r2.x += n.G.L.x;
}
}
- if (fabs(r1.y - r2.y) > n.G.L.y / 2) {
- if (r1.y < n.G.L.y / 2) {
- r2.y -= n.G.L.y;
+ if (n.G.edges[i].crossings[1]) {
+ if (r1.y < r2.y) {
+ r1.y += n.G.L.y;
} else {
r2.y += n.G.L.y;
}
@@ -111,20 +111,20 @@ void animate::post_fracture(network &n) {
graph::coordinate r1 = n.G.vertices[n.G.edges[i].v[0]].r;
graph::coordinate r2 = n.G.vertices[n.G.edges[i].v[1]].r;
- if (fabs(r1.x - r2.x) > n.G.L.x / 2) {
- if (r1.x < n.G.L.x / 2) {
- r1.x += n.G.L.x;
- } else {
- r2.x += n.G.L.x;
- }
- }
- if (fabs(r1.y - r2.y) > n.G.L.y / 2) {
- if (r1.y < n.G.L.y / 2) {
- r1.y += n.G.L.y;
- } else {
- r2.y += n.G.L.y;
- }
- }
+ if (n.G.edges[i].crossings[0]) {
+ if (r1.x < r2.x) {
+ r1.x += n.G.L.x;
+ } else {
+ r2.x += n.G.L.x;
+ }
+ }
+ if (n.G.edges[i].crossings[1]) {
+ if (r1.y < r2.y) {
+ r1.y += n.G.L.y;
+ } else {
+ r2.y += n.G.L.y;
+ }
+ }
glVertex2d(r1.x, r1.y);
glVertex2d(r2.x, r2.y);