many changes to introduce two-component, elastic-like fracture

1 files changed, 141 insertions, 63 deletions

diff --git a/lib/src/network.cpp b/lib/src/network.cpp
index 1edc973..8b0e8df 100644
--- a/lib/src/network.cpp
+++ b/lib/src/network.cpp
@@ -17,14 +17,19 @@ class nofuseException: public std::exception
   }
 } nofuseex;
 
-network::network(const graph& G, cholmod_common *c) : c(c), G(G), fuses(G.edges.size(), false), thresholds(G.edges.size(), 1) {
+problem::problem(const graph& G, unsigned axis, cholmod_sparse *vcmat, cholmod_common *c) : G(G), voltcurmat(vcmat), c(c), axis(axis) {
   b = CHOL_F(zeros)(G.vertices.size(), 1, CHOLMOD_REAL, c);
   for (unsigned i = 0; i < G.edges.size(); i++) {
-    double v0y = G.vertices[G.edges[i].v[0]].r.y;
-    double v1y = G.vertices[G.edges[i].v[1]].r.y;
-
-    if (G.edges[i].crossings[1]) {
-      bool ind = v0y < v1y ? 0 : 1;
+    graph::coordinate v0 = G.vertices[G.edges[i].v[0]].r;
+    graph::coordinate v1 = G.vertices[G.edges[i].v[1]].r;
+
+    if (G.edges[i].crossings[axis]) {
+      bool ind;
+      if (axis == 1) {
+        ind = v0.y < v1.y ? 0 : 1;
+      } else {
+        ind = v0.x < v1.x ? 0 : 1;
+      }
 
       ((double *)b->x)[G.edges[i].v[ind]] += 1.0;
       ((double *)b->x)[G.edges[i].v[!ind]] -= 1.0;
@@ -78,8 +83,10 @@ network::network(const graph& G, cholmod_common *c) : c(c), G(G), fuses(G.edges.
   factor = CHOL_F(analyze)(laplacian, c);
   CHOL_F(factorize)(laplacian, factor, c);
   CHOL_F(free_sparse)(&laplacian, c);
+}
 
-  t = CHOL_F(allocate_triplet)(G.edges.size(), G.vertices.size(), 2 * G.edges.size(), 0, CHOLMOD_REAL, c);
+problem::problem(const graph& G, unsigned axis, cholmod_common *c) : problem(G, axis, NULL, c) {
+  cholmod_triplet *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();
 
@@ -98,39 +105,70 @@ network::network(const graph& G, cholmod_common *c) : c(c), G(G), fuses(G.edges.
   CHOL_F(free_triplet)(&t, c);
 }
 
-network::network(const network &other) : c(other.c), G(other.G), fuses(other.fuses), thresholds(other.thresholds) {
+problem::problem(const problem& other) : G(other.G), c(other.c), axis(other.axis) {
   b = CHOL_F(copy_dense)(other.b, c);
   factor = CHOL_F(copy_factor)(other.factor, c);
   voltcurmat = CHOL_F(copy_sparse)(other.voltcurmat, c);
 } 
 
-network::~network() {
-  CHOL_F(free_sparse)(&voltcurmat, c);
+problem::~problem() {
   CHOL_F(free_dense)(&b, c);
   CHOL_F(free_factor)(&factor, c);
+  CHOL_F(free_sparse)(&voltcurmat, c);
 }
 
-void network::set_thresholds(double beta, std::mt19937& rng) {
-  if (beta == 0.0) {
-    /* zero beta doesn't make any sense computationally, we interpret it as "infinity" */
-    for (long double& threshold : thresholds) {
-      threshold = 1.0;
-    }
-  } else {
-    std::uniform_real_distribution<long double> d(0.0, 1.0);
+current_info problem::solve(std::vector<bool>& fuses) {
+  cholmod_dense *x = CHOL_F(solve)(CHOLMOD_A, factor, b, c);
 
-    for (long double& threshold : thresholds) {
-      threshold = std::numeric_limits<long double>::lowest();
+  if (((double *)x->x)[0] != ((double *)x->x)[0]) {
+    throw nanex;
+  }
 
-      while (threshold == std::numeric_limits<long double>::lowest()) {
-        threshold = logl(d(rng)) / (long double)beta;
+  cholmod_dense *y = CHOL_F(allocate_dense)(G.edges.size(), 1, G.edges.size(), CHOLMOD_REAL, c);
+
+  double alpha[2] = {1, 0};
+  double beta[2] = {0, 0};
+  CHOL_F(sdmult)(voltcurmat, 0, alpha, beta, x, y, c);
+
+  std::vector<double> currents(G.edges.size());
+
+  double total_current = 0;
+
+  for (int i = 0; i < G.edges.size(); i++) {
+    if (fuses[i]) {
+      currents[i] = 0;
+    } else {
+      currents[i] = ((double *)y->x)[i];
+
+      graph::coordinate v0 = G.vertices[G.edges[i].v[0]].r;
+      graph::coordinate v1 = G.vertices[G.edges[i].v[1]].r;
+
+      if (G.edges[i].crossings[axis]) {
+        bool comp;
+        if (axis == 1) {
+          comp = v0.y > v1.y;
+        } else {
+          comp = v0.x > v1.x;
+        }
+
+        if (comp) {
+          currents[i] += 1.0;
+          total_current += currents[i];
+        } else {
+          currents[i] -= 1.0;
+          total_current -= currents[i];
+        }
       }
     }
   }
+
+  CHOL_F(free_dense)(&x, c);
+  CHOL_F(free_dense)(&y, c);
+
+  return {total_current, currents};
 }
 
-void network::break_edge(unsigned e, bool unbreak) {
-  fuses[e] = !unbreak;
+void problem::break_edge(unsigned e, bool unbreak) {
   unsigned v0 = G.edges[e].v[0];
   unsigned v1 = G.edges[e].v[1];
 
@@ -167,78 +205,116 @@ void network::break_edge(unsigned e, bool unbreak) {
   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;
+  graph::coordinate r0 = G.vertices[v0].r;
+  graph::coordinate r1 = G.vertices[v1].r;
 
-  if (G.edges[e].crossings[1]) {
-    bool ind = v0y < v1y ? unbreak : !unbreak;
+  if (G.edges[e].crossings[axis]) {
+    bool ind;
+    if (axis == 1) {
+      ind = r0.y < r1.y ? unbreak : !unbreak;
+    } else {
+      ind = r0.x < r1.x ? unbreak : !unbreak;
+    }
 
     ((double *)b->x)[G.edges[e].v[ind]] -= 1.0;
     ((double *)b->x)[G.edges[e].v[!ind]] += 1.0;
   }
 }
 
-current_info network::get_current_info() {
-  cholmod_dense *x = CHOL_F(solve)(CHOLMOD_A, factor, b, c);
+network::network(const graph& G) : G(G), fuses(G.edges.size()), thresholds(G.edges.size()) {}
 
-  if (((double *)x->x)[0] != ((double *)x->x)[0]) {
-    throw nanex;
-  }
+network::network(const network& n) : G(n.G), fuses(n.fuses), thresholds(n.thresholds) {}
 
-  cholmod_dense *y = CHOL_F(allocate_dense)(G.edges.size(), 1, G.edges.size(), CHOLMOD_REAL, c);
+void network::set_thresholds(double beta, std::mt19937& rng) {
+  if (beta == 0.0) {
+    /* zero beta doesn't make any sense computationally, we interpret it as "infinity" */
+    for (long double& threshold : thresholds) {
+      threshold = 1.0;
+    }
+  } else {
+    std::uniform_real_distribution<long double> d(0.0, 1.0);
 
-  double alpha[2] = {1, 0};
-  double beta[2] = {0, 0};
-  CHOL_F(sdmult)(voltcurmat, 0, alpha, beta, x, y, c);
+    for (long double& threshold : thresholds) {
+      threshold = std::numeric_limits<long double>::lowest();
 
-  std::vector<double> currents(G.edges.size());
+      while (threshold == std::numeric_limits<long double>::lowest()) {
+        threshold = logl(d(rng)) / (long double)beta;
+      }
+    }
+  }
+}
 
-  double total_current = 0;
+fuse_network::fuse_network(const graph& G, cholmod_common *c) : network(G), ohm(G, 1, c) {
+}
 
-  for (int i = 0; i < G.edges.size(); i++) {
-    if (fuses[i]) {
-      currents[i] = 0;
-    } else {
-      currents[i] = ((double *)y->x)[i];
+void fuse_network::fracture(hooks& m, double cutoff) {
+  m.pre_fracture(*this);
 
-      double v0y = G.vertices[G.edges[i].v[0]].r.y;
-      double v1y = G.vertices[G.edges[i].v[1]].r.y;
+  while (true) {
+    current_info ci = ohm.solve(fuses);
 
-      if (G.edges[i].crossings[1]) {
-        if (v0y > v1y) {
-          currents[i] += 1.0;
-          total_current += currents[i];
-        } else {
-          currents[i] -= 1.0;
-          total_current -= currents[i];
+    if (ci.conductivity < cutoff * G.vertices.size()) {
+      break;
+    }
+
+    unsigned max_pos = UINT_MAX;
+    long double max_val = std::numeric_limits<long double>::lowest();
+
+    for (unsigned i = 0; i < G.edges.size(); i++) {
+      if (!fuses[i] && fabs(ci.currents[i]) > cutoff) {
+        long double val = logl(fabs(ci.currents[i])) - thresholds[i];
+        if (val > max_val) {
+          max_val = val;
+          max_pos = i;
         }
       }
     }
 
+    if (max_pos == UINT_MAX)  {
+      throw nofuseex;
+    }
+
+    fuses[max_pos] = true;
+    ohm.break_edge(max_pos);
+
+    m.bond_broken(*this, ci, max_pos);
   }
 
-  CHOL_F(free_dense)(&x, c);
-  CHOL_F(free_dense)(&y, c);
+  m.post_fracture(*this);
+}
 
-  return {total_current, currents};
+elastic_network::elastic_network(const graph& G, cholmod_common *c) : network(G), hook_x(G, 1, c), hook_y(G, 0, c) {
+}
+
+elastic_network::elastic_network(const elastic_network& n) : network(n), hook_x(n.hook_x), hook_y(n.hook_y) {
 }
 
-void network::fracture(hooks& m, double cutoff) {
+void elastic_network::fracture(hooks& m, double weight, double cutoff) {
   m.pre_fracture(*this);
 
   while (true) {
-    current_info ci = this->get_current_info();
+    current_info cx = hook_x.solve(fuses);
+    current_info cy = hook_y.solve(fuses);
 
-    if (ci.conductivity < cutoff * G.vertices.size()) {
+    current_info ctot;
+
+    ctot.conductivity = sqrt(pow((1 - weight) * cx.conductivity, 2) + pow(weight * cy.conductivity, 2));
+
+    if (ctot.conductivity < cutoff * G.vertices.size()) {
       break;
     }
 
+    ctot.currents.resize(G.edges.size());
+    for (unsigned i = 0; i < G.edges.size(); i++) {
+      ctot.currents[i] = sqrt(pow((1 - weight) * cx.currents[i], 2) + pow(weight * cy.currents[i], 2));
+    }
+
     unsigned max_pos = UINT_MAX;
     long double max_val = std::numeric_limits<long double>::lowest();
 
     for (unsigned i = 0; i < G.edges.size(); i++) {
-      if (!fuses[i] && fabs(ci.currents[i]) > cutoff) {
-        long double val = logl(fabs(ci.currents[i])) - thresholds[i];
+      if (!fuses[i] && fabs(ctot.currents[i]) > cutoff) {
+        long double val = logl(fabs(ctot.currents[i])) - thresholds[i];
         if (val > max_val) {
           max_val = val;
           max_pos = i;
@@ -250,9 +326,11 @@ void network::fracture(hooks& m, double cutoff) {
       throw nofuseex;
     }
 
-    this->break_edge(max_pos);
+    fuses[max_pos] = true;
+    hook_x.break_edge(max_pos);
+    hook_y.break_edge(max_pos);
 
-    m.bond_broken(*this, ci, max_pos);
+    m.bond_broken(*this, ctot, max_pos);
   }
 
   m.post_fracture(*this);