diff options
| -rw-r--r-- | src/measurements.cpp | 47 | ||||
| -rw-r--r-- | src/measurements.hpp | 2 |
2 files changed, 37 insertions, 12 deletions
diff --git a/src/measurements.cpp b/src/measurements.cpp index 17fbdee..368902a 100644 --- a/src/measurements.cpp +++ b/src/measurements.cpp @@ -1,6 +1,7 @@ #include "measurements.hpp" #include <iostream> +#include <cstdio> void update_distribution_file(std::string id, const std::vector<uint64_t>& data, unsigned 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"; @@ -64,15 +65,10 @@ void update_field_file(std::string id, const std::vector<T>& data, unsigned N, d file_out.close(); } -template <class T> -std::vector<fftw_complex> data_transform(unsigned Mx, unsigned My, const std::vector<T>& data, fftw_plan forward_plan, double *fftw_forward_in, fftw_complex *fftw_forward_out) { - for (unsigned i = 0; i < Mx * My; i++) { - fftw_forward_in[i] = (double)data[i]; - } - +fftw_complex* data_transform(unsigned Mx, unsigned My, fftw_plan forward_plan, double *fftw_forward_in, fftw_complex *fftw_forward_out) { fftw_execute(forward_plan); - std::vector<fftw_complex> output(Mx * (My / 2 + 1)); + fftw_complex* output = (fftw_complex*)malloc(Mx * (My / 2 + 1) * sizeof(fftw_complex)); for (unsigned i = 0; i < Mx * (My / 2 + 1); i++) { output[i][0] = fftw_forward_out[i][0]; @@ -83,7 +79,7 @@ std::vector<fftw_complex> data_transform(unsigned Mx, unsigned My, const std::ve } template <class T> -void correlation(unsigned Mx, unsigned 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) { +void correlation(unsigned Mx, unsigned My, std::vector<std::vector<T>>& data, const fftw_complex* tx1, const fftw_complex* tx2, fftw_plan reverse_plan, fftw_complex *fftw_reverse_in, double *fftw_reverse_out) { for (unsigned 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]; @@ -91,8 +87,10 @@ void correlation(unsigned Mx, unsigned My, std::vector<T>& data, const std::vect fftw_execute(reverse_plan); - for (unsigned 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)); + for (unsigned j = 0; j < data.size(); j++) { + for (unsigned i = 0; i < (Mx / 2 + 1) * (My / 2 + 1); i++) { + data[j][i] += (T)pow(fftw_reverse_out[Mx * (i / (Mx / 2 + 1)) + i % (Mx / 2 + 1)] / (Mx * My), j + 1); + } } } @@ -131,7 +129,9 @@ ma::ma(double Lx, double Ly, unsigned Mx, unsigned My, double beta, unsigned Ncu Caa(Ncum), Cdd(Ncum), Cll(Ncum), - Cee(Ncum) + Cee(Ncum), + CDD(Ncum), + CsD(Ncum) { N = 0; Nc = 0; @@ -146,6 +146,8 @@ ma::ma(double Lx, double Ly, unsigned Mx, unsigned My, double beta, unsigned Ncu Cdd[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0); Cll[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0); Cee[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0); + CDD[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0); + CsD[i].resize((Mx / 2 + 1) * (My / 2 + 1), 0); } // FFTW setup for correlation functions @@ -185,6 +187,8 @@ ma::~ma() { update_field_file("Caa", Caa, Na, Lx, Ly, beta, Mx, My); update_field_file("Cll", Cll, N, Lx, Ly, beta, Mx, My); update_field_file("Cee", Cee, N, Lx, Ly, beta, Mx, My); + update_field_file("CDD", CDD, N, Lx, Ly, beta, Mx, My); + update_field_file("CsD", CsD, N, Lx, Ly, beta, Mx, My); //stress_file.close(); } @@ -234,7 +238,9 @@ void ma::post_fracture(network &n) { fftw_forward_in[edge_r_to_ind(n.G.dual_vertices[n.G.dual_edges[edge].v[1]].r, Lx, Ly, Mx, My)] = 0.5; } - autocorrelation(Mx, My, Css, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); + fftw_complex *tss = data_transform(Mx, My, forward_plan, fftw_forward_in, fftw_forward_out); + + correlation(Mx, My, Css, tss, tss, reverse_plan, fftw_reverse_in, fftw_reverse_out); unsigned crack_component = component[n.G.dual_edges[crack.front()].v[0]]; @@ -289,6 +295,23 @@ void ma::post_fracture(network &n) { autocorrelation(Mx, My, Cee, forward_plan, fftw_forward_in, fftw_forward_out, reverse_plan, fftw_reverse_in, fftw_reverse_out); + /// damage at end + std::fill_n(fftw_forward_in, Mx * My, 0.0); + + for (unsigned i = 0; i < n.G.edges.size(); i++) { + if (n.fuses[i]) { + fftw_forward_in[edge_r_to_ind(n.G.edges[i].r, Lx, Ly, Mx, My)] += 1.0; + } + } + + fftw_complex *tDD = data_transform(Mx, My, forward_plan, fftw_forward_in, fftw_forward_out); + + correlation(Mx, My, CDD, tDD, tDD, reverse_plan, fftw_reverse_in, fftw_reverse_out); + correlation(Mx, My, CsD, tss, tDD, reverse_plan, fftw_reverse_in, fftw_reverse_out); + + free(tss); + free(tDD); + std::fill_n(fftw_forward_in, Mx * My, 0.0); // rewind the last avalanche diff --git a/src/measurements.hpp b/src/measurements.hpp index ba52720..a8251dd 100644 --- a/src/measurements.hpp +++ b/src/measurements.hpp @@ -48,6 +48,8 @@ class ma : public hooks { std::vector<std::vector<uint64_t>> Cdd; // damage-damage distribution std::vector<std::vector<uint64_t>> Cll; // damage-damage distribution std::vector<std::vector<uint64_t>> Cee; // damage-damage distribution + std::vector<std::vector<uint64_t>> CDD; // damage-damage distribution + std::vector<std::vector<uint64_t>> CsD; // damage-damage distribution uint64_t Nc; uint64_t NC; uint64_t Na; |