simplified random graph creation to only use uniform voronoi sites

4 files changed, 9 insertions, 82 deletions

diff --git a/lib/include/graph.hpp b/lib/include/graph.hpp
index ea08214..38c8edd 100644
--- a/lib/include/graph.hpp
+++ b/lib/include/graph.hpp
@@ -38,6 +38,6 @@ class graph {
     std::vector<edge> dual_edges;
 
     graph(unsigned int Nx, unsigned int Ny);
-    graph(double Lx, double Ly, std::mt19937& rng, double relax = 0.01, double step = 1.9);
+    graph(double Lx, double Ly, std::mt19937& rng);
 };
 
diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp
index 7c50ed0..2bd62a5 100644
--- a/lib/src/graph.cpp
+++ b/lib/src/graph.cpp
@@ -136,7 +136,7 @@ unsigned int get_triangle_signature(unsigned int j1, unsigned int j2, unsigned i
   }
 }
 
-graph::graph(double Lx, double Ly, std::mt19937& rng, double relax, double step) {
+graph::graph(double Lx, double Ly, std::mt19937& rng) {
   L = {Lx, Ly};
 
   // randomly choose N to be floor(Lx * Ly / 2) or ceil(Lx * Ly / 2) with
@@ -150,8 +150,8 @@ graph::graph(double Lx, double Ly, std::mt19937& rng, double relax, double step)
   // the coordinates of the lattice, from which a delaunay triangulation
   // 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)}};
+  for (vertex &v : vertices) {
+    v = {{Lx * d(rng), Ly * d(rng)}};
   }
 
   // set up the voronoi objects
@@ -160,72 +160,6 @@ graph::graph(double Lx, double Ly, std::mt19937& rng, double relax, double step)
   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
-    // the original points for a total of nine. note that the index of each
-    // point quotient 9 is equal to the original index (we'll use this to
-    // translate back)
-    for (unsigned int i = 0; i < nv; i++) {
-      const vertex& v = vertices[i];
-
-      points[9 * i + 0] = {v.r.x - L.x, v.r.y - L.y};
-      points[9 * i + 1] = {v.r.x + 0.0, v.r.y - L.y};
-      points[9 * i + 2] = {v.r.x + L.x, v.r.y - L.y};
-      points[9 * i + 3] = {v.r.x - L.x, v.r.y + 0.0};
-      points[9 * i + 4] = {v.r.x + 0.0, v.r.y + 0.0};
-      points[9 * i + 5] = {v.r.x + L.x, v.r.y + 0.0};
-      points[9 * i + 6] = {v.r.x - L.x, v.r.y + L.y};
-      points[9 * i + 7] = {v.r.x + 0.0, v.r.y + L.y};
-      points[9 * i + 8] = {v.r.x + L.x, v.r.y + L.y};
-    }
-
-    // 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 == 4) {
-        double Cx = 0.0;
-        double Cy = 0.0;
-        double A = 0.0;
-
-        const jcv_graphedge* e = site->edges;
-        while (e) {
-          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;
-        }
-
-        A /= 2;
-        Cx /= 6 * A;
-        Cy /= 6 * A;
-
-        double dx = Cx - vertices[ind / 9].r.x;
-        double dy = Cy - vertices[ind / 9].r.y;
-
-        double dist = pow(dx, 2) + pow(dy, 2);
-        if (dist > cur_diff) {
-          cur_diff = dist;
-        }
-        vertices[ind / 9] = {{mod(vertices[ind / 9].r.x + rstep * dx, L.x), mod(vertices[ind / 9].r.y + rstep * dy, L.y)}};
-      }
-    }
-    max_difference = cur_diff;
-    jcv_diagram_free(&diagram);
-    memset(&diagram, 0, sizeof(jcv_diagram));
-  }
-
   for (unsigned int i = 0; i < nv; i++) {
     const vertex& v = vertices[i];
     points[9 * i + 0] = {v.r.x - L.x, v.r.y - L.y};
@@ -238,6 +172,7 @@ graph::graph(double Lx, double Ly, std::mt19937& rng, double relax, double step)
     points[9 * i + 7] = {v.r.x + 0.0, 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(), &bounds, &diagram);
 
   const jcv_site* sites = jcv_diagram_get_sites(&diagram);
diff --git a/src/animate_fracture.cpp b/src/animate_fracture.cpp
index 62c5b1e..9ae79cc 100644
--- a/src/animate_fracture.cpp
+++ b/src/animate_fracture.cpp
@@ -34,9 +34,8 @@ int main(int argc, char* argv[]) {
   double Lx = 16.0;
   double Ly = 16.0;
   double  beta = 0.5;
-  double w = 0.01;
 
-  while ((opt = getopt(argc, argv, "X:Y:N:b:w:")) != -1) {
+  while ((opt = getopt(argc, argv, "X:Y:N:b:")) != -1) {
     switch (opt) {
       case 'N':
         N = (unsigned int)atof(optarg);
@@ -50,9 +49,6 @@ int main(int argc, char* argv[]) {
       case 'b':
         beta = atof(optarg);
         break;
-      case 'w':
-        w = atof(optarg);
-        break;
       default:
         exit(1);
     }
@@ -67,7 +63,7 @@ int main(int argc, char* argv[]) {
   std::mt19937 rng{seeds};
 
   for (unsigned int trial = 0; trial < N; trial++) {
-    graph G(Lx, Ly, rng, w);
+    graph G(Lx, Ly, rng);
     network network(G, &c);
     network.set_thresholds(beta, rng);
     network.fracture(meas);
diff --git a/src/fracture.cpp b/src/fracture.cpp
index 47b4856..3505166 100644
--- a/src/fracture.cpp
+++ b/src/fracture.cpp
@@ -34,9 +34,8 @@ 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:w:")) != -1) {
+  while ((opt = getopt(argc, argv, "N:X:Y:b:")) != -1) {
     switch (opt) {
       case 'N':
         N = (unsigned int)atof(optarg);
@@ -50,9 +49,6 @@ int main(int argc, char* argv[]) {
       case 'b':
         beta = atof(optarg);
         break;
-      case 'w':
-        w = atof(optarg);
-        break;
       default:
         exit(1);
     }
@@ -69,7 +65,7 @@ int main(int argc, char* argv[]) {
   for (unsigned int trial = 0; trial < N; trial++) {
     while (true) {
       try {
-        graph G(Lx, Ly, rng, w);
+        graph G(Lx, Ly, rng);
         network network(G, &c);
         network.set_thresholds(beta, rng);
         network.fracture(meas);