summaryrefslogtreecommitdiff
path: root/lib
diff options
context:
space:
mode:
Diffstat (limited to 'lib')
-rw-r--r--lib/include/current_info.hpp10
-rw-r--r--lib/include/hooks.hpp4
-rw-r--r--lib/include/network.hpp8
-rw-r--r--lib/src/network.cpp82
4 files changed, 32 insertions, 72 deletions
diff --git a/lib/include/current_info.hpp b/lib/include/current_info.hpp
new file mode 100644
index 0000000..90dd3c7
--- /dev/null
+++ b/lib/include/current_info.hpp
@@ -0,0 +1,10 @@
+
+#pragma once
+
+#include <vector>
+
+typedef struct current_info {
+ double conductivity;
+ std::vector<double> currents;
+} current_info;
+
diff --git a/lib/include/hooks.hpp b/lib/include/hooks.hpp
index 350f318..e617ac3 100644
--- a/lib/include/hooks.hpp
+++ b/lib/include/hooks.hpp
@@ -1,12 +1,14 @@
#pragma once
+#include "current_info.hpp"
+
class network;
class hooks {
public:
virtual void pre_fracture(const network&) {};
- virtual void bond_broken(const network&, const std::pair<double, std::vector<double>>&, unsigned int) {};
+ virtual void bond_broken(const network&, const current_info&, unsigned int) {};
virtual void post_fracture(network&) {}; // post fracture hook can be destructive
};
diff --git a/lib/include/network.hpp b/lib/include/network.hpp
index 5cbba0d..e6d67b9 100644
--- a/lib/include/network.hpp
+++ b/lib/include/network.hpp
@@ -6,10 +6,11 @@
#include <utility>
#include <random>
-#include "cholmod.h"
+#include <cholmod.h>
#include "graph.hpp"
#include "hooks.hpp"
+#include "current_info.hpp"
#ifdef FRACTURE_LONGINT
@@ -40,9 +41,8 @@ class network {
~network();
void set_thresholds(double, std::mt19937&);
- void break_edge(unsigned int);
- void add_edge(unsigned int);
- std::pair<double, std::vector<double>> currents();
+ void break_edge(unsigned int, bool unbreak = false);
+ current_info get_current_info();
void fracture(hooks&, double cutoff = 1e-14);
};
diff --git a/lib/src/network.cpp b/lib/src/network.cpp
index a27ca39..9f6cb0d 100644
--- a/lib/src/network.cpp
+++ b/lib/src/network.cpp
@@ -56,7 +56,6 @@ network::network(const graph& G, cholmod_common *c) : c(c), G(G), fuses(G.edges.
t = CHOL_F(allocate_triplet)(G.edges.size(), G.vertices.size(), 2 * G.edges.size(), 0, CHOLMOD_REAL, c);
-
t->nnz = 2 * G.edges.size();
for (unsigned int i = 0; i < G.edges.size(); i++) {
@@ -89,68 +88,17 @@ network::~network() {
void network::set_thresholds(double beta, std::mt19937& rng) {
std::uniform_real_distribution<long double> d(0.0, 1.0);
- for (unsigned int i = 0; i < G.edges.size(); i++) {
- long double x = 0.0;
+ for (long double& threshold : thresholds) {
+ threshold = 0.0;
- while (x == 0.0) {
- x = expl(logl(d(rng)) / (long double)beta);
+ while (threshold == 0.0) {
+ threshold = expl(logl(d(rng)) / (long double)beta);
}
-
- thresholds[i] = x;
- }
-}
-
-void network::add_edge(unsigned int e) {
- fuses[e] = false;
- unsigned int v0 = G.edges[e][0];
- unsigned int v1 = G.edges[e][1];
-
- unsigned int n = factor->n;
-
- cholmod_sparse *update_mat = CHOL_F(allocate_sparse)(n, n, 2, true, true, 0, CHOLMOD_REAL, c);
-
- unsigned int s1, s2;
- s1 = v0 < v1 ? v0 : v1;
- s2 = v0 < v1 ? v1 : v0;
-
- CHOL_INT *pp = (CHOL_INT *)update_mat->p;
- CHOL_INT *ii = (CHOL_INT *)update_mat->i;
- double *xx = (double *)update_mat->x;
-
- for (unsigned int i = 0; i <= s1; i++) {
- pp[i] = 0;
- }
-
- for (unsigned int i = s1 + 1; i <= n; i++) {
- pp[i] = 2;
- }
-
- ii[0] = s1;
- ii[1] = s2;
- xx[0] = -1;
- xx[1] = 1;
-
- cholmod_sparse *perm_update_mat = CHOL_F(submatrix)(
- update_mat, (CHOL_INT *)factor->Perm, factor->n, NULL, -1, true, true, c);
-
- CHOL_F(updown)(false, perm_update_mat, factor, c);
-
- CHOL_F(free_sparse)(&perm_update_mat, c);
- CHOL_F(free_sparse)(&update_mat, c);
-
- double v0y = G.vertices[v0].r.y;
- double v1y = G.vertices[v1].r.y;
-
- if (fabs(v0y - v1y) > 0.5) {
- bool ind = v0y < v1y ? 0 : 1;
-
- ((double *)b->x)[G.edges[e][ind]] += 1.0;
- ((double *)b->x)[G.edges[e][!ind]] -= 1.0;
}
}
-void network::break_edge(unsigned int e) {
- fuses[e] = true;
+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];
@@ -182,7 +130,7 @@ void network::break_edge(unsigned int e) {
cholmod_sparse *perm_update_mat = CHOL_F(submatrix)(
update_mat, (CHOL_INT *)factor->Perm, factor->n, NULL, -1, true, true, c);
- CHOL_F(updown)(false, perm_update_mat, factor, c);
+ CHOL_F(updown)(unbreak, perm_update_mat, factor, c);
CHOL_F(free_sparse)(&perm_update_mat, c);
CHOL_F(free_sparse)(&update_mat, c);
@@ -191,14 +139,14 @@ void network::break_edge(unsigned int e) {
double v1y = G.vertices[v1].r.y;
if (fabs(v0y - v1y) > 0.5) {
- bool ind = v0y < v1y ? 0 : 1;
+ bool ind = v0y < v1y ? unbreak : !unbreak;
((double *)b->x)[G.edges[e][ind]] -= 1.0;
((double *)b->x)[G.edges[e][!ind]] += 1.0;
}
}
-std::pair<double, std::vector<double>> network::currents() {
+current_info network::get_current_info() {
cholmod_dense *x = CHOL_F(solve)(CHOLMOD_A, factor, b, c);
if (((double *)x->x)[0] != ((double *)x->x)[0]) {
@@ -240,16 +188,16 @@ std::pair<double, std::vector<double>> network::currents() {
CHOL_F(free_dense)(&x, c);
CHOL_F(free_dense)(&y, c);
- return std::make_pair(total_current, currents);
+ return {total_current, currents};
}
void network::fracture(hooks& m, double cutoff) {
m.pre_fracture(*this);
while (true) {
- auto currents = this->currents();
+ current_info ci = this->get_current_info();
- if (currents.first < cutoff * G.vertices.size()) {
+ if (ci.conductivity < cutoff * G.vertices.size()) {
break;
}
@@ -257,8 +205,8 @@ void network::fracture(hooks& m, double cutoff) {
long double max_val = 0;
for (unsigned int i = 0; i < G.edges.size(); i++) {
- if (!fuses[i] && fabs(currents.second[i]) > cutoff) {
- long double val = (long double)fabs(currents.second[i]) / thresholds[i];
+ if (!fuses[i] && fabs(ci.currents[i]) > cutoff) {
+ long double val = (long double)fabs(ci.currents[i]) / thresholds[i];
if (val > max_val) {
max_val = val;
max_pos = i;
@@ -272,7 +220,7 @@ void network::fracture(hooks& m, double cutoff) {
this->break_edge(max_pos);
- m.bond_broken(*this, currents, max_pos);
+ m.bond_broken(*this, ci, max_pos);
}
m.post_fracture(*this);