diff options
| -rw-r--r-- | lib/src/network.cpp | 28 | ||||
| -rw-r--r-- | src/animate.cpp | 7 | ||||
| -rw-r--r-- | src/fracture.cpp | 21 | ||||
| -rw-r--r-- | src/measurements.cpp | 70 | ||||
| -rw-r--r-- | src/measurements.hpp | 4 | 
5 files changed, 54 insertions, 76 deletions
diff --git a/lib/src/network.cpp b/lib/src/network.cpp index dcc5e0e..ed23d45 100644 --- a/lib/src/network.cpp +++ b/lib/src/network.cpp @@ -90,6 +90,9 @@ void network::fracture(hooks& m, bool one_axis) {  }  void network::update_backbone(const std::vector<double>& c) { + +        bool done_x = px.sol.conductivity[0] < 1.0 / G.edges.size(); +        bool done_y = py.sol.conductivity[1] < 1.0 / G.edges.size();    class get_cycle_edges : public boost::default_dfs_visitor {      public:        std::set<unsigned> tE; @@ -191,7 +194,7 @@ void network::update_backbone(const std::vector<double>& c) {            other_cycles.push_back(tcycle);          } -        if (std::any_of(other_cycles.begin(), other_cycles.end(), [](std::array<unsigned, 2> c){return c[0] % 2 == 0 && c[1] % 2 == 0;})) { +        if (std::any_of(other_cycles.begin(), other_cycles.end(), [done_x, done_y](std::array<unsigned, 2> c){return (c[0] % 2 == 0 && c[1] % 2 == 0) || ((c[0] % 2 == 0 && done_x) || (c[1] % 2 == 0 && done_y));})) {            seen_pairs[{G.edges[i].v[0], G.edges[i].v[1]}] = true;            backbone[i] = true;            boost::remove_edge(G.edges[i].v[0], G.edges[i].v[1], bG); @@ -304,6 +307,8 @@ void network::update_backbone(const std::vector<double>& c) {            bool any_intervening_e_1 = false;            bool any_intervening2 = false;            bool any_intervening_e_2 = false; +          unsigned ie1 = 0; +          unsigned ie2 = 0;            for (unsigned k = il + 1; k < ig; k++) {              if (!boost::same_component(G.vertices[v].nd[i], G.vertices[v].nd[k], ds)) { @@ -320,33 +325,24 @@ void network::update_backbone(const std::vector<double>& c) {            if (any_intervening2 && !any_intervening1) {              for (unsigned k = il + 1; k <= ig; k++) {                if (!backbone[G.vertices[v].ne[k]]) { -                any_intervening_e_1 = true; -                break; +                ie1++;                }              }            }            if (any_intervening1 && !any_intervening2) {              for (unsigned k = ig + 1; k%G.vertices[v].nd.size() != il + 1; k++) {                if (!backbone[G.vertices[v].ne[k%G.vertices[v].nd.size()]]) { -                any_intervening_e_2 = true; -                break; +                ie2++;                }              }            } -          if ((any_intervening1 && any_intervening2) || (any_intervening1 && any_intervening_e_2) || (any_intervening2 && any_intervening_e_1)) { +          if ((any_intervening1 && any_intervening2) || (any_intervening1 && ie2 > 1) || (any_intervening2 && ie1 > 1)) {              found_pair = true;              p1 = il;              p2 = ig;              d1 = G.vertices[v].nd[il];              d2 = G.vertices[v].nd[ig]; -            break; -          } -        } -        } -      } -    } -    if (found_pair) {        std::list<unsigned> cedges = {};        get_cycle_edges vis(cedges);        std::vector<boost::default_color_type> cm1(G.dual_vertices.size()); @@ -390,7 +386,7 @@ void network::update_backbone(const std::vector<double>& c) {            other_cycles.push_back(tcycle);          } -        if (std::any_of(other_cycles.begin(), other_cycles.end(), [](std::array<unsigned, 2> c){return c[0] % 2 == 0 && c[1] % 2 == 0;})) { +        if (std::any_of(other_cycles.begin(), other_cycles.end(), [done_x, done_y](std::array<unsigned, 2> c){return ((c[0] % 2 == 0 && c[1] % 2 == 0) || (c[0] % 2 == 0 && done_x)) || (c[1] % 2 == 0 && done_y);})) {            seen_pairs[{d1, d2}] = true;        std::set<unsigned> left_edges = {}; @@ -456,6 +452,10 @@ void network::update_backbone(const std::vector<double>& c) {        }    }    } +        } +      } +    } +  }  }  void network::break_edge(unsigned e, bool unbreak) { diff --git a/src/animate.cpp b/src/animate.cpp index 532bd70..b1031dc 100644 --- a/src/animate.cpp +++ b/src/animate.cpp @@ -89,12 +89,14 @@ void animate::bond_broken(const network& n, const current_info& cur, unsigned i)      */     bool weird = false; +   unsigned nw = 0;      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; +    nw++;      draw = true;        } else if (n.backbone[j] && !n.fuses[j] && j != i) {      glColor3f(0.8f, 0.8f, 0.8f); // Blue @@ -102,6 +104,9 @@ void animate::bond_broken(const network& n, const current_info& cur, unsigned i)        } else if (!n.fuses[j]) {      glColor3f(0.0f, 0.0f, 0.0f); // Blue      draw = true; +      } else if (i == j) { +    glColor3f(0.0f, 1.0f, 0.0f); // Blue +    draw = true;        }        if (draw) { @@ -130,7 +135,7 @@ void animate::bond_broken(const network& n, const current_info& cur, unsigned i)        }   glEnd();   glFlush(); - if (weird) {std::cout << "\n"; getchar();} + if (nw > 2) {std::cout << "\n"; getchar();}  }  void animate::post_fracture(network &n) { diff --git a/src/fracture.cpp b/src/fracture.cpp index 483a3d2..5ba60e2 100644 --- a/src/fracture.cpp +++ b/src/fracture.cpp @@ -38,9 +38,10 @@ int main(int argc, char* argv[]) {    unsigned n = 128;    double a = 1.0;    bool use_aN = false; +  bool one_side = false;    double w = 0.5; -  while ((opt = getopt(argc, argv, "N:X:Y:b:n:a:w:")) != -1) { +  while ((opt = getopt(argc, argv, "N:X:Y:b:n:a:w:o")) != -1) {      switch (opt) {        case 'N':          N = (unsigned)atof(optarg); @@ -65,6 +66,9 @@ int main(int argc, char* argv[]) {        case 'w':          w = atof(optarg);          break; +      case 'o': +        one_side = true; +        break;        default:          exit(1);      } @@ -77,10 +81,7 @@ int main(int argc, char* argv[]) {    std::mt19937 rng{seeds};    if (use_aN) { -    unsigned Mx = (unsigned)4*sqrt(2*n*a); -    unsigned My = (unsigned)4*sqrt(2*n/a); - -    ma meas(n, a, beta); +    ma meas(n, a, beta, w, one_side);      for (unsigned trial = 0; trial < N; trial++) {        while (true) { @@ -88,7 +89,7 @@ int main(int argc, char* argv[]) {            graph G(n, a, rng);            elastic_network fuse_network(G, &c, w);            fuse_network.set_thresholds(beta, rng); -          fuse_network.fracture(meas, false); +          fuse_network.fracture(meas, one_side);            break;          } catch (std::exception &e) {            std::cout << e.what() << '\n'; @@ -99,11 +100,7 @@ int main(int argc, char* argv[]) {          break;      }    } else { -    // fourier transforms of powers of two are faster -    unsigned Mx = pow(2, ceil(log2(4*Lx))); -    unsigned My = pow(2, ceil(log2(4*Ly))); - -    ma meas(Lx, Ly, beta); +    ma meas(Lx, Ly, beta, w, one_side);      for (unsigned trial = 0; trial < N; trial++) {        while (true) { @@ -111,7 +108,7 @@ int main(int argc, char* argv[]) {            graph G(Lx, Ly);            elastic_network fuse_network(G, &c, w);            fuse_network.set_thresholds(beta, rng); -          fuse_network.fracture(meas, false); +          fuse_network.fracture(meas, one_side);            break;          } catch (std::exception &e) {            std::cout << e.what() << '\n'; diff --git a/src/measurements.cpp b/src/measurements.cpp index c58fd84..96b4ccb 100644 --- a/src/measurements.cpp +++ b/src/measurements.cpp @@ -142,45 +142,55 @@ unsigned edge_r_to_ind(graph::coordinate r, double Lx, double Ly, unsigned Mx, u    return floor((Mx * r.x) / Lx) + Mx * floor((My * r.y) / Ly);  } -ma::ma(unsigned n, double a, double beta) : +ma::ma(unsigned n, double a, double beta, double weight, bool one) :    G(2 * n),    sn(2 * n),    ss(2 * n),    sm(2 * n),    sl(2 * n), -  sb(n), +  sb(n + 1),    sd(3 * n), -  sc(2 * n),    sa(3 * n),    sA(3 * n),    si(10000),    sI(10000)  {    if (beta != 0.0) { -    model_string = "fracture_" + std::to_string(n) + "_" + std::to_string(a) + "_" + std::to_string(beta) + "_"; +    model_string = "fracture_" + std::to_string(n) + "_" + std::to_string(a) + "_" + std::to_string(beta) + "_" + std::to_string(weight) + "_";    } else { -    model_string = "fracture_" + std::to_string(n) + "_" + std::to_string(a) + "_INF_"; +    model_string = "fracture_" + std::to_string(n) + "_" + std::to_string(a) + "_INF_" + std::to_string(weight) + "_"; +  } + +  if (one) { +    model_string = model_string + "o_"; +  } else { +    model_string = model_string + "t_";    }  } -ma::ma(unsigned Lx, unsigned Ly, double beta) : +ma::ma(unsigned Lx, unsigned Ly, double beta, double weight, bool one) :    G(Lx * Ly / 2),    sn(Lx * Ly / 2),    ss(Lx * Ly / 2),    sm(Lx * Ly / 2),    sl(Lx * Ly / 2), -  sb(Lx * Ly / 2), +  sb(Lx * Ly / 2 + 1),    sd(Lx * Ly), -  sc(Lx * Ly / 2),    sa(Lx * Ly),    sA(Lx * Ly),    si(10000),    sI(10000)  {    if (beta != 0.0) { -    model_string = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_"; +    model_string = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_" + std::to_string(weight) + "_";    } else { -    model_string = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_INF_"; +    model_string = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_INF_" + std::to_string(weight) + "_"; +  } + +  if (one) { +    model_string = model_string + "o_"; +  } else { +    model_string = model_string + "t_";    }  } @@ -191,7 +201,6 @@ ma::~ma() {    update_distribution_file("sl", sl, model_string);    update_distribution_file("sb", sb, model_string);    update_distribution_file("sd", sd, model_string); -  update_distribution_file("sc", sc, model_string);    update_distribution_file("sa", sa, model_string);    update_distribution_file("sA", sA, model_string);    update_distribution_file("si", si, model_string); @@ -260,43 +269,10 @@ void ma::post_fracture(network &n) {      sn[components[i].size() - 1]++;    } -  auto av_it = avalanches.rbegin(); - -  while (true) { -    for (unsigned e : *av_it) { -      boost::remove_edge(n.G.dual_edges[e].v[0], n.G.dual_edges[e].v[1], G); -      n.break_edge(e, true); -    } - -    auto cracks = find_minimal_crack(G, n); - -    if (cracks.size() == 0) { -      break; -    } - -    av_it++; -  } - -  num = boost::connected_components(G, &component[0]); - -  std::vector<std::list<graph::coordinate>> new_components(num); - -  for (unsigned i = 0; i < n.G.dual_vertices.size(); i++) { -    new_components[component[i]].push_back(n.G.dual_vertices[i].r); -  } - -  for (unsigned i = 0; i < num; i++) { -    sc[new_components[i].size() - 1]++; -  } - -  /* -  current_info ct = n.get_current_info(); - -    std::vector<bool> vertex_in(n.G.vertices.size());    for (unsigned i = 0; i < n.G.edges.size(); i++) { -    if (ct.currents[i] > CURRENT_CUTOFF) { +    if (!n.backbone[i]) {        vertex_in[n.G.edges[i].v[0]] = true;        vertex_in[n.G.edges[i].v[1]] = true;      } @@ -308,9 +284,9 @@ void ma::post_fracture(network &n) {      if (vertex_in[i]) bb_size++;    } -  sb[bb_size - 1]++; -  */ +  sb[bb_size]++; +  auto av_it = avalanches.rbegin();    av_it++;    while (av_it != avalanches.rend()) { diff --git a/src/measurements.hpp b/src/measurements.hpp index 274a550..c211974 100644 --- a/src/measurements.hpp +++ b/src/measurements.hpp @@ -44,8 +44,8 @@ class ma : public hooks {      std::list<unsigned> last_avalanche;      std::string model_string; -    ma(unsigned Lx, unsigned Ly, double beta); -    ma(unsigned n, double a, double beta); +    ma(unsigned Lx, unsigned Ly, double beta, double w, bool o); +    ma(unsigned n, double a, double beta, double w, bool o);      ~ma();      void pre_fracture(const network &) override;  | 
