cleaned up graph a bit, redefined assert so that certain run errors in jc_voronoi can be caught

4 files changed, 54 insertions, 18 deletions

diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp
index 3d5647f..b9d4f03 100644
--- a/lib/src/graph.cpp
+++ b/lib/src/graph.cpp
@@ -6,8 +6,12 @@
 #define JCV_REAL_TYPE double
 #define JCV_ATAN2 atan2
 #define JCV_FLT_MAX 1.7976931348623157E+308
+#undef assert
+#define assert_error_report_helper(cond) "assertion failed: "
+#define assert(cond) {if(!(cond)) { std::cerr << assert_error_report_helper(cond) "\n"; throw assert_error_report_helper(cond); } }
 #include <jc_voronoi.h>
 
+// actual mod for floats
 double mod(double a, double b) {
   if (a >= 0) {
     return fmod(a, b);
@@ -95,26 +99,30 @@ class eulerException: public std::exception
 graph::graph(double Lx, double Ly, std::mt19937& rng, double relax, double step) {
   L = {Lx, Ly};
 
+  // randomly choose N to be floor(Lx * Ly / 2) or ceil(Lx * Ly / 2) with
+  // probability proportional to the distance from each
   std::uniform_real_distribution<double> d(0.0, 1.0);
   unsigned int N = round(Lx * Ly / 2 + d(rng) - 0.5);
 
   unsigned int nv = N;
-
   vertices.resize(nv);
 
   // the coordinates of the lattice, from which a delaunay triangulation
-  // and corresponding voronoi tessellation will be built
+  // and corresponding voronoi tessellation will be built. Everyone is in the
+  // rectangle (0, 0) (Lx, Ly)
   for (unsigned int i = 0; i < nv; i++) {
     vertices[i] = {{L.x * d(rng), L.y * d(rng)}};
   }
 
-  double max_difference = std::numeric_limits<double>::max();
+  // set up the voronoi objects
   jcv_diagram diagram;
   memset(&diagram, 0, sizeof(jcv_diagram));
+  jcv_rect bounds = {{-Lx, -Ly}, {2 * Lx, 2 * Ly}};
   std::vector<jcv_point> points(9 * nv);
 
   double rstep = sqrt(step);
 
+  double max_difference = std::numeric_limits<double>::max();
   while (max_difference > pow(relax, 2) * 2 * N) {
     double cur_diff = 0;
     // to make the resulting tessellation periodic, we tile eight (!) copies of
@@ -123,6 +131,7 @@ graph::graph(double Lx, double Ly, std::mt19937& rng, double relax, double step)
     // translate back)
     for (unsigned int i = 0; i < nv; i++) {
       const vertex& v = vertices[i];
+
       points[9 * i + 0] = {v.r.x + 0.0, v.r.y + 0.0};
       points[9 * i + 1] = {v.r.x + L.x, v.r.y + 0.0};
       points[9 * i + 2] = {v.r.x - L.x, v.r.y + 0.0};
@@ -134,27 +143,33 @@ graph::graph(double Lx, double Ly, std::mt19937& rng, double relax, double step)
       points[9 * i + 8] = {v.r.x - L.x, v.r.y - L.y};
     }
 
-    jcv_diagram_generate(9 * nv, points.data(), NULL, &diagram);
+    // run voronoi
+    jcv_diagram_generate(9 * nv, points.data(), &bounds, &diagram);
+
+    // relax the sites by moving the vertices towards their centroids
     const jcv_site* sites = jcv_diagram_get_sites(&diagram);
-    
     for (int i = 0; i < diagram.numsites; i++) {
       const jcv_site* site = &sites[i];
       unsigned int ind = site->index;
       if (ind % 9 == 0) {
         double Cx = 0.0;
         double Cy = 0.0;
-        unsigned int ne = 0;
+        double A = 0.0;
 
         const jcv_graphedge* e = site->edges;
         while (e) {
-          ne++;
-          Cx += e->pos[0].x;
-          Cy += e->pos[0].y;
+          double dA = e->pos[0].x * e->pos[1].y - e->pos[1].x * e->pos[0].y;
+
+          A += dA;
+          Cx += (e->pos[0].x + e->pos[1].x) * dA;
+          Cy += (e->pos[0].y + e->pos[1].y) * dA;
+
           e = e->next;
         }
 
-        Cx /= ne;
-        Cy /= ne;
+        A /= 2;
+        Cx /= 6 * A;
+        Cy /= 6 * A;
 
         double dx = Cx - vertices[ind / 9].r.x;
         double dy = Cy - vertices[ind / 9].r.y;
@@ -183,7 +198,7 @@ graph::graph(double Lx, double Ly, std::mt19937& rng, double relax, double step)
     points[9 * i + 7] = {v.r.x + L.x, v.r.y - L.y};
     points[9 * i + 8] = {v.r.x - L.x, v.r.y - L.y};
   }
-  jcv_diagram_generate(9 * nv, points.data(), NULL, &diagram);
+  jcv_diagram_generate(9 * nv, points.data(), &bounds, &diagram);
 
   const jcv_site* sites = jcv_diagram_get_sites(&diagram);
 
diff --git a/src/fracture.cpp b/src/fracture.cpp
index 742a3f9..47b4856 100644
--- a/src/fracture.cpp
+++ b/src/fracture.cpp
@@ -34,8 +34,9 @@ int main(int argc, char* argv[]) {
   double Lx = 16;
   double Ly = 16;
   double  beta = 0.5;
+  double w = 0.01;
 
-  while ((opt = getopt(argc, argv, "N:X:Y:b:")) != -1) {
+  while ((opt = getopt(argc, argv, "N:X:Y:b:w:")) != -1) {
     switch (opt) {
       case 'N':
         N = (unsigned int)atof(optarg);
@@ -49,6 +50,9 @@ int main(int argc, char* argv[]) {
       case 'b':
         beta = atof(optarg);
         break;
+      case 'w':
+        w = atof(optarg);
+        break;
       default:
         exit(1);
     }
@@ -65,12 +69,13 @@ int main(int argc, char* argv[]) {
   for (unsigned int trial = 0; trial < N; trial++) {
     while (true) {
       try {
-        graph G(Lx, Ly, rng);
+        graph G(Lx, Ly, rng, w);
         network network(G, &c);
         network.set_thresholds(beta, rng);
         network.fracture(meas);
         break;
       } catch (std::exception &e) {
+        std::cout << e.what() << '\n';
       }
     }
 
diff --git a/src/measurements.cpp b/src/measurements.cpp
index 27f239d..15d37d1 100644
--- a/src/measurements.cpp
+++ b/src/measurements.cpp
@@ -112,6 +112,8 @@ ma::ma(double Lx, double Ly, unsigned int Mx, unsigned int My, double beta) :
   Lx(Lx), Ly(Ly), Mx(Mx), My(My), beta(beta), G(2 * (unsigned int)ceil(Lx * Ly / 2)),
   sc(2 * (unsigned int)ceil(Lx * Ly / 2), 0),
   sa(2 * (unsigned int)ceil(Lx * Ly / 2), 0),
+  sC(2 * (unsigned int)ceil(Lx * Ly / 2), 0),
+  sA(2 * (unsigned int)ceil(Lx * Ly / 2), 0),
   sd(2 * (unsigned int)ceil(Lx * Ly / 2), 0),
   sb(log2(Mx < My ? Mx : My) + 1, 0),
   Ccc((Mx / 2 + 1) * (My / 2 + 1), 0),
@@ -125,6 +127,8 @@ ma::ma(double Lx, double Ly, unsigned int Mx, unsigned int My, double beta) :
   N = 0;
   Nc = 0;
   Na = 0;
+  NC = 0;
+  NA = 0;
 
   // FFTW setup for correlation functions
   fftw_set_timelimit(1);
@@ -150,6 +154,8 @@ ma::~ma() {
 
   update_distribution_file("sa", sa, Na, Lx, Ly, beta);
   update_distribution_file("sc", sc, Nc, Lx, Ly, beta);
+  update_distribution_file("sA", sA, NA, Lx, Ly, beta);
+  update_distribution_file("sC", sC, NC, Lx, Ly, beta);
   update_distribution_file("sd", sd, N, Lx, Ly, beta);
   update_distribution_file("sb", sb, N, Lx, Ly, beta);
 
@@ -172,7 +178,9 @@ void ma::bond_broken(const network& net, const current_info& cur, unsigned int 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]++;
+    sA[avalanches.back().size() - 1]++;
     Na++;
+    NA++;
 
     std::fill_n(fftw_forward_in, Mx * My, 0.0);
 
@@ -224,8 +232,10 @@ void ma::post_fracture(network &n) {
       }
 
       sc[components[i].size() - 1]++;
+      sC[components[i].size() - 1]++;
       autocorrelation(Mx, My, Ccc, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
       Nc++;
+      NC++;
 
       for (auto it = components[i].begin(); it != components[i].end(); it++) {
         fftw_forward_in[edge_r_to_ind(n.G.dual_vertices[*it].r, Lx, Ly, Mx, My)] = 0.0;
@@ -252,10 +262,10 @@ void ma::post_fracture(network &n) {
   // spanning cluster
   std::fill_n(fftw_forward_in, Mx * My, 0.0);
 
-  for (unsigned int i = 0; i < n.G.dual_vertices.size(); i++) {
-    if (component[i] == crack_component)  {
-      fftw_forward_in[edge_r_to_ind(n.G.dual_vertices[i].r, Lx, Ly, Mx, My)] += 1.0;
-    }
+  sC[components[crack_component].size() - 1]++;
+  NC++;
+  for (auto it = components[crack_component].begin(); it != components[crack_component].end(); it++) {
+    fftw_forward_in[edge_r_to_ind(n.G.dual_vertices[*it].r, Lx, Ly, Mx, My)] += 1.0;
   }
 
   autocorrelation(Mx, My, Cmm, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out);
@@ -277,6 +287,8 @@ void ma::post_fracture(network &n) {
   std::fill_n(fftw_forward_in, Mx * My, 0.0);
 
   // rewind the last avalanche
+  sA[avalanches.back().size() - 1]++;
+  NA++;
   for (auto e : avalanches.back()) {
     boost::remove_edge(n.G.dual_edges[e].v[0], n.G.dual_edges[e].v[1], G);
     n.break_edge(e, true);
diff --git a/src/measurements.hpp b/src/measurements.hpp
index bd93cfc..3d4ff3e 100644
--- a/src/measurements.hpp
+++ b/src/measurements.hpp
@@ -36,6 +36,8 @@ class ma : public hooks {
     // measurement storage
     std::vector<uint64_t> sc; // cluster size distribution
     std::vector<uint64_t> sa; // avalanche size distribution
+    std::vector<uint64_t> sC; // cluster size distribution
+    std::vector<uint64_t> sA; // avalanche size distribution
     std::vector<uint64_t> sd; // avalanche size distribution
     std::vector<uint64_t> sb; // bin size distribution
     std::vector<uint64_t> Ccc; // cluster-cluster correlations
@@ -46,7 +48,9 @@ class ma : public hooks {
     std::vector<uint64_t> Cll; // damage-damage distribution
     std::vector<uint64_t> Cee; // damage-damage distribution
     uint64_t Nc;
+    uint64_t NC;
     uint64_t Na;
+    uint64_t NA;
 
   public:
     long double lv;