added field for crack backbone/damage correlations

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;