mostly implemented the ability to find dead bonds using topological properties

7 files changed, 94 insertions, 40 deletions

diff --git a/src/analysis_tools.hpp b/src/analysis_tools.hpp
index 659cacf..c3bffa2 100644
--- a/src/analysis_tools.hpp
+++ b/src/analysis_tools.hpp
@@ -1,10 +1,4 @@
 
-#include <boost/graph/adjacency_list.hpp>
-#include <boost/graph/connected_components.hpp>
-#include <boost/graph/depth_first_search.hpp>
-#include <boost/range/combine.hpp>
-#include <boost/foreach.hpp>
-
 #include <vector>
 #include <algorithm>
 #include <cmath>
@@ -13,14 +7,6 @@
 
 #include <network.hpp>
 
-struct EdgeProperties {
-  unsigned index;
-};
-
-typedef boost::adjacency_list<boost::listS, boost::vecS, boost::undirectedS, boost::no_property, EdgeProperties> Graph;
-typedef boost::graph_traits<Graph>::vertex_descriptor Vertex;
-typedef boost::graph_traits<Graph>::vertices_size_type VertexIndex;
-
 
 template<class T>
 bool is_shorter(const std::list<T> &, const std::list<T> &);
diff --git a/src/animate.cpp b/src/animate.cpp
index 9ba335b..532bd70 100644
--- a/src/animate.cpp
+++ b/src/animate.cpp
@@ -19,9 +19,10 @@ void animate::pre_fracture(const network &n) {
   boost::remove_edge_if(trivial, G);
 
   glClearColor(1.0f, 1.0f, 1.0f, 1.0f );
+  glLineWidth(5);
   glClear(GL_COLOR_BUFFER_BIT);
   glBegin(GL_LINES);
-  glColor3f(0.0f, 0.0f, 0.0f);
+  glColor3f(0.9f, 0.9f, 0.9f);
   for (unsigned i = 0; i < n.G.edges.size(); i++) {
     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;
@@ -59,19 +60,62 @@ void animate::bond_broken(const network& n, const current_info& cur, unsigned i)
 
   boost::add_edge(n.G.dual_edges[i].v[0], n.G.dual_edges[i].v[1], {i}, G);
 
+  glClearColor(1.0f, 1.0f, 1.0f, 1.0f );
+  glClear(GL_COLOR_BUFFER_BIT);
    glBegin(GL_LINES);          // Each set of 3 vertices form a triangle
-    glColor3f(1.0f, 1.0f, 1.0f); // Blue
-    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 (n.G.edges[i].crossings[0]) {
+  /*
+    glColor3f(0.0f, 0.0f, 0.0f); // Blue
+    graph::coordinate r1 = n.G.dual_vertices[n.G.dual_edges[i].v[0]].r;
+    graph::coordinate r2 = n.G.dual_vertices[n.G.dual_edges[i].v[1]].r;
+
+    if (n.G.dual_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 (n.G.dual_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);
+    */
+
+   bool weird = false;
+
+    for (unsigned j = 0; j < n.G.edges.size(); j++) {
+      bool draw = false;
+      if (cur.currents[j] < 1e-9 && !n.backbone[j]) {
+    glColor3f(1.0f, 0.0f, 0.0f); // Blue
+    weird = true;
+    draw = true;
+      } else if (n.backbone[j] && !n.fuses[j] && j != i) {
+    glColor3f(0.8f, 0.8f, 0.8f); // Blue
+    draw = true;
+      } else if (!n.fuses[j]) {
+    glColor3f(0.0f, 0.0f, 0.0f); // Blue
+    draw = true;
+      }
+
+      if (draw) {
+    graph::coordinate r1 = n.G.vertices[n.G.edges[j].v[0]].r;
+    graph::coordinate r2 = n.G.vertices[n.G.edges[j].v[1]].r;
+
+    if (n.G.edges[j].crossings[0]) {
+      if (r1.x < r2.x) {
+        r1.x += n.G.L.x;
+      } else {
+        r2.x += n.G.L.x;
+      }
+    }
+    if (n.G.edges[j].crossings[1]) {
       if (r1.y < r2.y) {
         r1.y += n.G.L.y;
       } else {
@@ -81,12 +125,17 @@ void animate::bond_broken(const network& n, const current_info& cur, unsigned i)
 
     glVertex2d(r1.x, r1.y);
     glVertex2d(r2.x, r2.y);
+      }
+
+      }
  glEnd();
  glFlush();
+ if (weird) {std::cout << "\n"; getchar();}
 }
 
 void animate::post_fracture(network &n) {
 
+  /*
   std::list<unsigned> crack;
 //  unsigned crack_component = component[n.G.dual_edges[crack.front()].v[0]];
   unsigned crack_component = 10000;
@@ -223,5 +272,6 @@ void animate::post_fracture(network &n) {
         glFlush();
     }
   }
+*/
 }
 
diff --git a/src/animate_fracture.cpp b/src/animate_fracture.cpp
index c3b4a69..9d875e6 100644
--- a/src/animate_fracture.cpp
+++ b/src/animate_fracture.cpp
@@ -31,20 +31,20 @@ int main(int argc, char* argv[]) {
   int opt;
 
   unsigned N = 1;
-  double Lx = 16.0;
-  double Ly = 16.0;
+  unsigned n = 16;
+  double a = 16.0;
   double  beta = 0.5;
 
-  while ((opt = getopt(argc, argv, "X:Y:N:b:")) != -1) {
+  while ((opt = getopt(argc, argv, "n:a:N:b:")) != -1) {
     switch (opt) {
       case 'N':
         N = (unsigned)atof(optarg);
         break;
-      case 'X':
-        Lx = atof(optarg);
+      case 'n':
+        n = atoi(optarg);
         break;
-      case 'Y':
-        Ly = atof(optarg);
+      case 'a':
+        a = atof(optarg);
         break;
       case 'b':
         beta = atof(optarg);
@@ -57,16 +57,16 @@ int main(int argc, char* argv[]) {
   cholmod_common c;
   CHOL_F(start)(&c);
 
-  animate meas(Lx, Ly, 700, argc, argv);
+  animate meas(sqrt(2*n *a), sqrt( 2*n / a), 700, argc, argv);
 
   randutils::auto_seed_128 seeds;
   std::mt19937 rng{seeds};
 
   for (unsigned trial = 0; trial < N; trial++) {
-    graph G(Lx, Ly, rng);
+    graph G(n, a, rng);
     elastic_network elastic_network(G, &c);
     elastic_network.set_thresholds(beta, rng);
-    elastic_network.fracture(meas);
+    elastic_network.fracture(meas, true);
 
     if (quit.load())
       break;
diff --git a/src/animate_fracture_square.cpp b/src/animate_fracture_square.cpp
index 84cc8a5..bc4c3a2 100644
--- a/src/animate_fracture_square.cpp
+++ b/src/animate_fracture_square.cpp
@@ -68,7 +68,7 @@ int main(int argc, char* argv[]) {
   for (unsigned trial = 0; trial < N; trial++) {
     elastic_network tmp_network(perm_network);
     tmp_network.set_thresholds(beta, rng);
-    tmp_network.fracture(meas);
+    tmp_network.fracture(meas, false);
 
     if (quit.load())
       break;
diff --git a/src/fracture.cpp b/src/fracture.cpp
index 7df3f67..483a3d2 100644
--- a/src/fracture.cpp
+++ b/src/fracture.cpp
@@ -38,8 +38,9 @@ int main(int argc, char* argv[]) {
   unsigned n = 128;
   double a = 1.0;
   bool use_aN = false;
+  double w = 0.5;
 
-  while ((opt = getopt(argc, argv, "N:X:Y:b:n:a:")) != -1) {
+  while ((opt = getopt(argc, argv, "N:X:Y:b:n:a:w:")) != -1) {
     switch (opt) {
       case 'N':
         N = (unsigned)atof(optarg);
@@ -61,6 +62,9 @@ int main(int argc, char* argv[]) {
         a = atof(optarg);
         use_aN = true;
         break;
+      case 'w':
+        w = atof(optarg);
+        break;
       default:
         exit(1);
     }
@@ -82,9 +86,9 @@ int main(int argc, char* argv[]) {
       while (true) {
         try {
           graph G(n, a, rng);
-          elastic_network fuse_network(G, &c);
+          elastic_network fuse_network(G, &c, w);
           fuse_network.set_thresholds(beta, rng);
-          fuse_network.fracture(meas);
+          fuse_network.fracture(meas, false);
           break;
         } catch (std::exception &e) {
           std::cout << e.what() << '\n';
@@ -105,9 +109,9 @@ int main(int argc, char* argv[]) {
       while (true) {
         try {
           graph G(Lx, Ly);
-          elastic_network fuse_network(G, &c);
+          elastic_network fuse_network(G, &c, w);
           fuse_network.set_thresholds(beta, rng);
-          fuse_network.fracture(meas);
+          fuse_network.fracture(meas, false);
           break;
         } catch (std::exception &e) {
           std::cout << e.what() << '\n';
diff --git a/src/measurements.cpp b/src/measurements.cpp
index 0732d24..c58fd84 100644
--- a/src/measurements.cpp
+++ b/src/measurements.cpp
@@ -153,7 +153,8 @@ ma::ma(unsigned n, double a, double beta) :
   sc(2 * n),
   sa(3 * n),
   sA(3 * n),
-  si(10000)
+  si(10000),
+  sI(10000)
 {
   if (beta != 0.0) {
     model_string = "fracture_" + std::to_string(n) + "_" + std::to_string(a) + "_" + std::to_string(beta) + "_";
@@ -173,7 +174,8 @@ ma::ma(unsigned Lx, unsigned Ly, double beta) :
   sc(Lx * Ly / 2),
   sa(Lx * Ly),
   sA(Lx * Ly),
-  si(10000)
+  si(10000),
+  sI(10000)
 {
   if (beta != 0.0) {
     model_string = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_";
@@ -193,6 +195,7 @@ ma::~ma() {
   update_distribution_file("sa", sa, model_string);
   update_distribution_file("sA", sA, model_string);
   update_distribution_file("si", si, model_string);
+  update_distribution_file("sI", sI, model_string);
 }
 
 void ma::pre_fracture(const network&) {
@@ -212,8 +215,10 @@ void ma::bond_broken(const network& net, const current_info& cur, unsigned i) {
   }
 
   for (unsigned j = 0; j < cur.currents.size(); j++) {
-    if (logl(cur.currents[j]) >= -100 && logl(cur.currents[j]) < 0) {
+    if (logl(cur.currents[j]) >= -100 && logl(cur.currents[j]) < 0 && (!net.backbone[j] || j == i)) {
       si[(unsigned)(10000 * (logl(cur.currents[j]) + 100) / 100)]++;
+    } else if (logl(cur.currents[j]) >= -100 && logl(cur.currents[j]) < 0 && net.backbone[j] && j != i) {
+      sI[(unsigned)(10000 * (logl(cur.currents[j]) + 100) / 100)]++;
     }
   }
 
@@ -226,7 +231,13 @@ void ma::post_fracture(network &n) {
   std::vector<unsigned> component(boost::num_vertices(G));
   unsigned num = boost::connected_components(G, &component[0]);
   if (post_cracks.size() > 2 || post_cracks.size() == 0) {
-    std::cout << post_cracks.size() << "\n";
+    for (auto c : post_cracks) {
+      for (unsigned e : c.second) {
+        std::cout << e << " ";
+      }
+      std::cout << "\n";
+    }
+      getchar();
     throw badcycleex;
   }
   for (auto c : post_cracks) {
@@ -278,6 +289,7 @@ void ma::post_fracture(network &n) {
     sc[new_components[i].size() - 1]++;
   }
 
+  /*
   current_info ct = n.get_current_info();
 
 
@@ -297,6 +309,7 @@ void ma::post_fracture(network &n) {
   }
 
   sb[bb_size - 1]++;
+  */
 
   av_it++;
 
diff --git a/src/measurements.hpp b/src/measurements.hpp
index 5b76e26..274a550 100644
--- a/src/measurements.hpp
+++ b/src/measurements.hpp
@@ -35,6 +35,7 @@ class ma : public hooks {
     std::vector<uint64_t> sA; // non-final avalanche size distribution
 
     std::vector<uint64_t> si;
+    std::vector<uint64_t> sI;
 
   public:
     long double lv;