diff options
| -rw-r--r-- | lib/include/network.hpp | 1 | ||||
| -rw-r--r-- | lib/src/network.cpp | 10 | ||||
| -rw-r--r-- | src/measurements.cpp | 6 |
3 files changed, 9 insertions, 8 deletions
diff --git a/lib/include/network.hpp b/lib/include/network.hpp index 849975f..d955e43 100644 --- a/lib/include/network.hpp +++ b/lib/include/network.hpp @@ -5,6 +5,7 @@ #include <functional> #include <utility> #include <random> +#include <limits> #include <cholmod.h> diff --git a/lib/src/network.cpp b/lib/src/network.cpp index bc7a0c6..9cb1007 100644 --- a/lib/src/network.cpp +++ b/lib/src/network.cpp @@ -89,10 +89,10 @@ void network::set_thresholds(double beta, std::mt19937& rng) { std::uniform_real_distribution<long double> d(0.0, 1.0); for (long double& threshold : thresholds) { - threshold = 0.0; + threshold = std::numeric_limits<long double>::lowest(); - while (threshold == 0.0) { - threshold = expl(logl(d(rng)) / (long double)beta); + while (threshold == std::numeric_limits<long double>::lowest()) { + threshold = logl(d(rng)) / (long double)beta; } } } @@ -202,11 +202,11 @@ void network::fracture(hooks& m, double cutoff) { } unsigned int max_pos = UINT_MAX; - long double max_val = 0; + long double max_val = std::numeric_limits<long double>::lowest(); for (unsigned int i = 0; i < G.edges.size(); i++) { if (!fuses[i] && fabs(ci.currents[i]) > cutoff) { - long double val = (long double)fabs(ci.currents[i]) / thresholds[i]; + long double val = logl(fabs(ci.currents[i])) - thresholds[i]; if (val > max_val) { max_val = val; max_pos = i; diff --git a/src/measurements.cpp b/src/measurements.cpp index 5dc6a64..98795f1 100644 --- a/src/measurements.cpp +++ b/src/measurements.cpp @@ -103,11 +103,11 @@ 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.vertices[n.G.dual_edges[edge][0]].r.x - n.G.vertices[n.G.dual_edges[edge][1]].r.x); + double dx = fabs(n.G.dual_vertices[n.G.dual_edges[edge][0]].r.x - n.G.dual_vertices[n.G.dual_edges[edge][1]].r.x); if (dx > n.G.L.x / 2) { crossing_count[0]++; } - double dy = fabs(n.G.vertices[n.G.dual_edges[edge][0]].r.y - n.G.vertices[n.G.dual_edges[edge][1]].r.y); + double dy = fabs(n.G.dual_vertices[n.G.dual_edges[edge][0]].r.y - n.G.dual_vertices[n.G.dual_edges[edge][1]].r.y); if (dy > n.G.L.y / 2) { crossing_count[1]++; } @@ -292,7 +292,7 @@ void ma::pre_fracture(const network &) { } void ma::bond_broken(const network& net, const current_info& cur, unsigned int i) { - long double c = cur.conductivity / fabs(cur.currents[i]) * net.thresholds[i]; + long double c = logl(cur.conductivity / fabs(cur.currents[i])) + net.thresholds[i]; if (c > lv && avalanches.back().size() > 0) { sa[avalanches.back().size() - 1]++; Na++; |