summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--src/measurements.cpp47
-rw-r--r--src/measurements.hpp2
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;