From 9e1610143f0e96b77ca962a7113d79a2fbcbf6f5 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Fri, 28 Dec 2018 17:33:35 -0500 Subject: partially fixed problems that arise in small systems by reworking the way that boundary edges are identified --- lib/include/array_hash.hpp | 27 +++++++++++++++++++++++++++ lib/include/graph.hpp | 2 ++ lib/include/network.hpp | 1 + lib/src/graph.cpp | 34 ++++++++++------------------------ lib/src/network.cpp | 33 +++++++++++++++++++++------------ 5 files changed, 61 insertions(+), 36 deletions(-) create mode 100644 lib/include/array_hash.hpp (limited to 'lib') 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 +#include + +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; + } + }; +} + 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 #include +#include "array_hash.hpp" class graph { public: @@ -25,6 +26,7 @@ class graph { typedef struct edge { std::array v; coordinate r; + std::array 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 // actual mod for floats @@ -65,27 +64,6 @@ graph::graph(unsigned int Nx, unsigned int Ny) { } } -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() @@ -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, 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]; -- cgit v1.2.3-54-g00ecf