summaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/animate.cpp64
-rw-r--r--src/animate.hpp4
-rw-r--r--src/animate_fracture.cpp16
-rw-r--r--src/fracture.cpp12
-rw-r--r--src/measurements.cpp41
-rw-r--r--src/measurements.hpp6
6 files changed, 73 insertions, 70 deletions
diff --git a/src/animate.cpp b/src/animate.cpp
index 8c2d12f..26f4996 100644
--- a/src/animate.cpp
+++ b/src/animate.cpp
@@ -1,15 +1,15 @@
#include "animate.hpp"
-animate::animate(unsigned int Lx, unsigned int Ly, unsigned int window_size, int argc, char *argv[]) : Lx(Lx), Ly(Ly), G(Lx * Ly) {
+animate::animate(double Lx, double Ly, unsigned int window_size, int argc, char *argv[]) : G(2 * (unsigned int)ceil(Lx * Ly / 2)) {
glutInit(&argc, argv);
glutInitDisplayMode(GLUT_SINGLE | GLUT_RGB);
- glutInitWindowSize(window_size, (unsigned int )(window_size * ((double)Ly / (double)Lx)));
+ glutInitWindowSize((unsigned int)(Lx / Ly * window_size), window_size);
glutCreateWindow("wolff");
glClearColor(0.0,0.0,0.0,0.0);
glMatrixMode(GL_PROJECTION);
glLoadIdentity();
- gluOrtho2D(-1.0, Lx + 3.0, -1.0 , Ly + 3.0);
+ gluOrtho2D(-1, Lx + 3, -1 , Ly + 3);
}
void animate::pre_fracture(const network &n) {
@@ -25,18 +25,18 @@ void animate::pre_fracture(const network &n) {
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 (fabs(r1.x - r2.x) > Lx / 2) {
- if (r1.x < Lx / 2) {
- r1.x += Lx;
+ if (fabs(r1.x - r2.x) > n.G.L.x / 2) {
+ if (r1.x < n.G.L.x / 2) {
+ r1.x += n.G.L.x;
} else {
- r2.x += Lx;
+ r2.x += n.G.L.x;
}
}
- if (fabs(r1.y - r2.y) > Ly / 2) {
- if (r1.y < Ly / 2) {
- r1.y += Ly;
+ if (fabs(r1.y - r2.y) > n.G.L.y / 2) {
+ if (r1.y < n.G.L.y / 2) {
+ r1.y += n.G.L.y;
} else {
- r2.y += Ly;
+ r2.y += n.G.L.y;
}
}
@@ -47,8 +47,8 @@ void animate::pre_fracture(const network &n) {
glFlush();
}
-void animate::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];
+void animate::bond_broken(const network& n, const current_info& cur, unsigned int i) {
+ long double c = logl(cur.conductivity / fabs(cur.currents[i])) + n.thresholds[i];
if (c > lv && avalanches.back().size() > 0) {
lv = c;
avalanches.push_back({i});
@@ -56,25 +56,25 @@ void animate::bond_broken(const network& net, const current_info& cur, unsigned
avalanches.back().push_back(i);
}
- boost::add_edge(net.G.dual_edges[i].v[0], net.G.dual_edges[i].v[1], {i}, G);
+ boost::add_edge(n.G.dual_edges[i].v[0], n.G.dual_edges[i].v[1], {i}, G);
glBegin(GL_LINES); // Each set of 3 vertices form a triangle
glColor3f(1.0f, 1.0f, 1.0f); // Blue
- graph::coordinate r1 = net.G.vertices[net.G.edges[i].v[0]].r;
- graph::coordinate r2 = net.G.vertices[net.G.edges[i].v[1]].r;
+ 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 (fabs(r1.x - r2.x) > Lx / 2) {
- if (r1.x < Lx / 2) {
- r2.x -= Lx;
+ if (fabs(r1.x - r2.x) > n.G.L.x / 2) {
+ if (r1.x < n.G.L.x / 2) {
+ r2.x -= n.G.L.x;
} else {
- r2.x += Lx;
+ r2.x += n.G.L.x;
}
}
- if (fabs(r1.y - r2.y) > Ly / 2) {
- if (r1.y < Ly / 2) {
- r2.y -= Ly;
+ if (fabs(r1.y - r2.y) > n.G.L.y / 2) {
+ if (r1.y < n.G.L.y / 2) {
+ r2.y -= n.G.L.y;
} else {
- r2.y += Ly;
+ r2.y += n.G.L.y;
}
}
@@ -111,18 +111,18 @@ void animate::post_fracture(network &n) {
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 (fabs(r1.x - r2.x) > Lx / 2) {
- if (r1.x < Lx / 2) {
- r1.x += Lx;
+ if (fabs(r1.x - r2.x) > n.G.L.x / 2) {
+ if (r1.x < n.G.L.x / 2) {
+ r1.x += n.G.L.x;
} else {
- r2.x += Lx;
+ r2.x += n.G.L.x;
}
}
- if (fabs(r1.y - r2.y) > Ly / 2) {
- if (r1.y < Ly / 2) {
- r1.y += Ly;
+ if (fabs(r1.y - r2.y) > n.G.L.y / 2) {
+ if (r1.y < n.G.L.y / 2) {
+ r1.y += n.G.L.y;
} else {
- r2.y += Ly;
+ r2.y += n.G.L.y;
}
}
diff --git a/src/animate.hpp b/src/animate.hpp
index 793771d..423eba7 100644
--- a/src/animate.hpp
+++ b/src/animate.hpp
@@ -9,14 +9,12 @@
class animate : public hooks {
private:
- unsigned int Lx;
- unsigned int Ly;
Graph G;
public:
long double lv;
std::list<std::list<unsigned int>> avalanches;
- animate(unsigned int Lx, unsigned int Ly, unsigned int window_size, int argc, char *argv[]);
+ animate(double Lx, double Ly, unsigned int window_size, int argc, char *argv[]);
void pre_fracture(const network &) override;
void bond_broken(const network& net, const current_info& cur, unsigned int i) override;
diff --git a/src/animate_fracture.cpp b/src/animate_fracture.cpp
index 10b63d0..62c5b1e 100644
--- a/src/animate_fracture.cpp
+++ b/src/animate_fracture.cpp
@@ -31,24 +31,28 @@ int main(int argc, char* argv[]) {
int opt;
unsigned int N = 1;
- unsigned int Lx = 16;
- unsigned int Ly = 16;
+ double Lx = 16.0;
+ double Ly = 16.0;
double beta = 0.5;
+ double w = 0.01;
- while ((opt = getopt(argc, argv, "N:X:Y:b:")) != -1) {
+ while ((opt = getopt(argc, argv, "X:Y:N:b:w:")) != -1) {
switch (opt) {
case 'N':
N = (unsigned int)atof(optarg);
break;
case 'X':
- Lx = atoi(optarg);
+ Lx = atof(optarg);
break;
case 'Y':
- Ly = atoi(optarg);
+ Ly = atof(optarg);
break;
case 'b':
beta = atof(optarg);
break;
+ case 'w':
+ w = atof(optarg);
+ break;
default:
exit(1);
}
@@ -63,7 +67,7 @@ int main(int argc, char* argv[]) {
std::mt19937 rng{seeds};
for (unsigned int trial = 0; trial < N; trial++) {
- graph G(Lx, Ly, rng, 0.4);
+ graph G(Lx, Ly, rng, w);
network network(G, &c);
network.set_thresholds(beta, rng);
network.fracture(meas);
diff --git a/src/fracture.cpp b/src/fracture.cpp
index 4b6fb5b..742a3f9 100644
--- a/src/fracture.cpp
+++ b/src/fracture.cpp
@@ -31,8 +31,8 @@ int main(int argc, char* argv[]) {
int opt;
unsigned int N = 1;
- unsigned int Lx = 16;
- unsigned int Ly = 16;
+ double Lx = 16;
+ double Ly = 16;
double beta = 0.5;
while ((opt = getopt(argc, argv, "N:X:Y:b:")) != -1) {
@@ -41,10 +41,10 @@ int main(int argc, char* argv[]) {
N = (unsigned int)atof(optarg);
break;
case 'X':
- Lx = atoi(optarg);
+ Lx = atof(optarg);
break;
case 'Y':
- Ly = atoi(optarg);
+ Ly = atof(optarg);
break;
case 'b':
beta = atof(optarg);
@@ -57,7 +57,7 @@ int main(int argc, char* argv[]) {
cholmod_common c;
CHOL_F(start)(&c);
- ma meas(Lx, Ly, Lx, Ly, beta);
+ ma meas(Lx, Ly, ceil(Lx), ceil(Ly), beta);
randutils::auto_seed_128 seeds;
std::mt19937 rng{seeds};
@@ -65,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, 0.4);
+ graph G(Lx, Ly, rng);
network network(G, &c);
network.set_thresholds(beta, rng);
network.fracture(meas);
diff --git a/src/measurements.cpp b/src/measurements.cpp
index 43483be..27f239d 100644
--- a/src/measurements.cpp
+++ b/src/measurements.cpp
@@ -2,7 +2,7 @@
#include "measurements.hpp"
#include <iostream>
-void update_distribution_file(std::string id, const std::vector<uint64_t>& data, unsigned int N, unsigned int Lx, unsigned int Ly, double beta) {
+void update_distribution_file(std::string id, const std::vector<uint64_t>& data, unsigned int N, double Lx, double Ly, double beta) {
std::string filename = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_" + id + ".dat";
std::ifstream file(filename);
@@ -31,7 +31,7 @@ void update_distribution_file(std::string id, const std::vector<uint64_t>& data,
}
template <class T>
-void update_field_file(std::string id, const std::vector<T>& data, unsigned int N, unsigned int Lx, unsigned int Ly, double beta, unsigned int Mx, unsigned int My) {
+void update_field_file(std::string id, const std::vector<T>& data, unsigned int N, double Lx, double Ly, double beta, unsigned int Mx, unsigned int My) {
std::string filename = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" + std::to_string(beta) + "_" + id + "_" + std::to_string(Mx) + "_" + std::to_string(My) + ".dat";
std::ifstream file(filename);
@@ -57,16 +57,16 @@ void update_field_file(std::string id, const std::vector<T>& data, unsigned int
}
template <class T>
-std::vector<fftw_complex> data_transform(unsigned int Lx, unsigned int Ly, const std::vector<T>& data, fftw_plan forward_plan, double *fftw_forward_in, fftw_complex *fftw_forward_out) {
- for (unsigned int i = 0; i < Lx * Ly; i++) {
+std::vector<fftw_complex> data_transform(unsigned int Mx, unsigned int My, const std::vector<T>& data, fftw_plan forward_plan, double *fftw_forward_in, fftw_complex *fftw_forward_out) {
+ for (unsigned int i = 0; i < Mx * My; i++) {
fftw_forward_in[i] = (double)data[i];
}
fftw_execute(forward_plan);
- std::vector<fftw_complex> output(Lx * (Ly / 2 + 1));
+ std::vector<fftw_complex> output(Mx * (My / 2 + 1));
- for (unsigned int i = 0; i < Lx * (Ly / 2 + 1); i++) {
+ for (unsigned int i = 0; i < Mx * (My / 2 + 1); i++) {
output[i][0] = fftw_forward_out[i][0];
output[i][1] = fftw_forward_out[i][1];
}
@@ -75,43 +75,44 @@ std::vector<fftw_complex> data_transform(unsigned int Lx, unsigned int Ly, const
}
template <class T>
-void correlation(unsigned int Lx, unsigned int Ly, std::vector<T>& data, const std::vector<fftw_complex>& tx1, const std::vector<fftw_complex>& tx2, fftw_plan reverse_plan, fftw_complex *fftw_reverse_in, double *fftw_reverse_out) {
- for (unsigned int i = 0; i < Lx * (Ly / 2 + 1); i++) {
+void correlation(unsigned int Mx, unsigned int My, std::vector<T>& data, const std::vector<fftw_complex>& tx1, const std::vector<fftw_complex>& tx2, fftw_plan reverse_plan, fftw_complex *fftw_reverse_in, double *fftw_reverse_out) {
+ for (unsigned int i = 0; i < Mx * (My / 2 + 1); i++) {
fftw_reverse_in[i][0] = tx1[i][0] * tx2[i][0] + tx1[i][1] * tx2[i][1];
fftw_reverse_in[i][1] = tx1[i][0] * tx2[i][1] - tx1[i][1] * tx2[i][0];
}
fftw_execute(reverse_plan);
- for (unsigned int i = 0; i < (Lx / 2 + 1) * (Ly / 2 + 1); i++) {
- data[i] += (T)(fftw_reverse_out[Lx * (i / (Lx / 2 + 1)) + i % (Lx / 2 + 1)] / (Lx * Ly));
+ for (unsigned int i = 0; i < (Mx / 2 + 1) * (My / 2 + 1); i++) {
+ data[i] += (T)(fftw_reverse_out[Mx * (i / (Mx / 2 + 1)) + i % (Mx / 2 + 1)] / (Mx * My));
}
}
template <class T>
-void autocorrelation(unsigned int Lx, unsigned int Ly, std::vector<T>& out_data, fftw_plan forward_plan, double *fftw_forward_in, fftw_complex *fftw_forward_out, fftw_plan reverse_plan, fftw_complex *fftw_reverse_in, double *fftw_reverse_out) {
+void autocorrelation(unsigned int Mx, unsigned int My, std::vector<T>& out_data, fftw_plan forward_plan, double *fftw_forward_in, fftw_complex *fftw_forward_out, fftw_plan reverse_plan, fftw_complex *fftw_reverse_in, double *fftw_reverse_out) {
fftw_execute(forward_plan);
- for (unsigned int i = 0; i < Ly * (Lx / 2 + 1); i++) {
+ for (unsigned int i = 0; i < My * (Mx / 2 + 1); i++) {
fftw_reverse_in[i][0] = pow(fftw_forward_out[i][0], 2) + pow(fftw_forward_out[i][1], 2);
fftw_reverse_in[i][1] = 0.0;
}
fftw_execute(reverse_plan);
- for (unsigned int i = 0; i < (Lx / 2 + 1) * (Ly / 2 + 1); i++) {
- out_data[i] += (T)(fftw_reverse_out[Lx * (i / (Lx / 2 + 1)) + i % (Lx / 2 + 1)] / (Lx * Ly));
+ for (unsigned int i = 0; i < (Mx / 2 + 1) * (My / 2 + 1); i++) {
+ out_data[i] += (T)(fftw_reverse_out[Mx * (i / (Mx / 2 + 1)) + i % (Mx / 2 + 1)] / (Mx * My));
}
}
-unsigned int edge_r_to_ind(graph::coordinate r, unsigned int Lx, unsigned int Ly, unsigned int Mx, unsigned int My) {
+unsigned int edge_r_to_ind(graph::coordinate r, double Lx, double Ly, unsigned int Mx, unsigned int My) {
return floor((Mx * r.x) / Lx) + Mx * floor((My * r.y) / Ly);
}
-ma::ma(unsigned int Lx, unsigned int Ly, unsigned int Mx, unsigned int My, double beta) : Lx(Lx), Ly(Ly), Mx(Mx), My(My), beta(beta), G(Lx * Ly),
- sc(Lx * Ly, 0),
- sa(Lx * Ly, 0),
- sd(Lx * Ly, 0),
+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),
+ 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),
Css((Mx / 2 + 1) * (My / 2 + 1), 0),
@@ -217,7 +218,7 @@ void ma::post_fracture(network &n) {
// non-spanning clusters
std::fill_n(fftw_forward_in, Mx * My, 0.0);
for (unsigned int i = 0; i < num; i++) {
- if (i != crack_component) {
+ if (i != crack_component && components[i].size() > 0) {
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)] += 1.0;
}
diff --git a/src/measurements.hpp b/src/measurements.hpp
index 685806b..bd93cfc 100644
--- a/src/measurements.hpp
+++ b/src/measurements.hpp
@@ -26,8 +26,8 @@ class ma : public hooks {
fftw_plan forward_plan;
fftw_plan reverse_plan;
unsigned int N;
- unsigned int Lx;
- unsigned int Ly;
+ double Lx;
+ double Ly;
unsigned int Mx;
unsigned int My;
double beta;
@@ -54,7 +54,7 @@ class ma : public hooks {
std::list<std::list<unsigned int>> avalanches;
- ma(unsigned int Lx, unsigned int Ly, unsigned int Mx, unsigned int My, double beta);
+ ma(double Lx, double Ly, unsigned int Mx, unsigned int My, double beta);
~ma();
void pre_fracture(const network &) override;