automated a field class

2 files changed, 57 insertions, 16 deletions

diff --git a/src/measurements.cpp b/src/measurements.cpp
index 15658d2..f42b5c7 100644
--- a/src/measurements.cpp
+++ b/src/measurements.cpp
@@ -62,6 +62,38 @@ void update_field_file(std::string id, uint64_t N, const std::vector<T>& data,
 
   file_out.close();
 }
+
+void update_field_file(std::string id, const RealField& data,
+                       std::string model_string) {
+  std::string filename = model_string + id + ".dat";
+  std::ifstream file(filename);
+
+  std::vector<long unsigned> data_old_N(data.N.size(), 0);
+  std::vector<double> data_old_F(data.F.size(), 0);
+
+  if (file.is_open()) {
+    for (unsigned j = 0; j < data.N.size(); j++) {
+      file >> data_old_N[j];
+    }
+    for (unsigned j = 0; j < data.F.size(); j++) {
+      file >> data_old_F[j];
+    }
+    file.close();
+  }
+
+  std::ofstream file_out(filename);
+
+  for (unsigned j = 0; j < data.N.size(); j++) {
+    file_out << data_old_N[j] + data.N[j] << " ";
+  }
+  file_out << "\n";
+  for (unsigned j = 0; j < data.F.size(); j++) {
+    file_out << data_old_F[j] + data.F[j] << " ";
+  }
+
+  file_out.close();
+}
+
 /*
 fftw_complex* data_transform(unsigned Mx, unsigned My, fftw_plan forward_plan,
 double *fftw_forward_in, fftw_complex *fftw_forward_out) {
@@ -171,7 +203,7 @@ ma::ma(unsigned n, double a, double beta, double weight)
       cm(pow(1 + 2 * (unsigned)ceil(sqrt(n)), 2)), cl(pow(1 + 2 * (unsigned)ceil(sqrt(n)), 2)),
       cb(pow(1 + 2 * (unsigned)ceil(sqrt(n)), 2)), ca(pow(1 + 2 * (unsigned)ceil(sqrt(n)), 2)),
       cA(pow(1 + 2 * (unsigned)ceil(sqrt(n)), 2)), cq(pow(1 + 2 * (unsigned)ceil(sqrt(n)), 2)),
-      cS(pow(4 * (unsigned)ceil(sqrt(n)), 2) / 2), NS(pow(4 * (unsigned)ceil(sqrt(n)), 2) / 2), last_clusters(2 * n) {
+      fS({sqrt(2 * n * a), sqrt(2 * n / a) / 2}, {4 * (unsigned)ceil(sqrt(n)), 2 * (unsigned)ceil(sqrt(n))}), last_clusters(2 * n) {
   if (beta != 0.0) {
     model_string = "fracture_" + std::to_string(n) + "_" + std::to_string(a) + "_" +
                    std::to_string(beta) + "_" + std::to_string(weight) + "_";
@@ -199,7 +231,8 @@ ma::ma(unsigned Lx, unsigned Ly, double beta, double weight)
       cs((Lx / 2 + 1) * (Ly / 2 + 1)), cm((Lx / 2 + 1) * (Ly / 2 + 1)),
       cl((Lx / 2 + 1) * (Ly / 2 + 1)), cb((Lx / 2 + 1) * (Ly / 2 + 1)),
       ca((Lx / 2 + 1) * (Ly / 2 + 1)), cA((Lx / 2 + 1) * (Ly / 2 + 1)),
-      cq((Lx / 2 + 1) * (Ly / 2 + 1)), cS(Lx * Ly / 2), NS(Lx * Ly / 2), last_clusters(Lx * Ly / 2) {
+      cq((Lx / 2 + 1) * (Ly / 2 + 1)), fS({(double)Lx, (double)Ly / 2}, {Lx, Ly
+          / 2}), last_clusters(Lx * Ly / 2) {
   if (beta != 0.0) {
     model_string = "fracture_" + std::to_string(Lx) + "_" + std::to_string(Ly) + "_" +
                    std::to_string(beta) + "_" + std::to_string(weight) + "_";
@@ -242,8 +275,7 @@ ma::~ma() {
   update_field_file("ca", Na, ca, model_string);
   update_field_file("cA", NA, cA, model_string);
   update_field_file("cq", Nq, cq, model_string);
-  update_field_file("cS", 0, cS, model_string);
-  update_field_file("NS", 0, NS, model_string);
+  update_field_file("fS", fS, model_string);
 }
 
 void ma::pre_fracture(const network&) {
@@ -514,16 +546,8 @@ void ma::post_fracture(network& n) {
           Δyr = Δy_tmpr > n.G.L.x / 2 ? n.G.L.x - Δy_tmpr : Δy_tmpr;
         }
 
-        cS[(unsigned)floor(Mx * (Δyl / n.G.L.y)) +
-           (Mx / 2) * (unsigned)floor(My * (Δxl / n.G.L.x))] +=
-            last_currents.currents[i] / last_currents.conductivity[crack.first[0] % 2];
-        NS[(unsigned)floor(Mx * (Δyl / n.G.L.y)) +
-           (Mx / 2) * (unsigned)floor(My * (Δxl / n.G.L.x))]++;
-        NS[(unsigned)floor(Mx * (Δyr / n.G.L.y)) +
-           (Mx / 2) * (unsigned)floor(My * (Δxr / n.G.L.x))]++;
-        cS[(unsigned)floor(Mx * (Δyr / n.G.L.y)) +
-           (Mx / 2) * (unsigned)floor(My * (Δxr / n.G.L.x))] +=
-            last_currents.currents[i] / last_currents.conductivity[crack.first[0] % 2];
+        fS.add({Δxr, Δyr}, last_currents.currents[i] / last_currents.conductivity[crack.first[0] % 2]);
+        fS.add({Δxl, Δyl}, last_currents.currents[i] / last_currents.conductivity[crack.first[0] % 2]);
       }
     }
 
diff --git a/src/measurements.hpp b/src/measurements.hpp
index 7bc8137..18b2e8e 100644
--- a/src/measurements.hpp
+++ b/src/measurements.hpp
@@ -12,6 +12,24 @@
 
 #include "analysis_tools.hpp"
 
+class RealField {
+  public:
+    coordinate L;
+    std::array<unsigned, 2> M;
+    std::vector<long unsigned> N;
+    std::vector<double> F;
+
+    RealField(coordinate L, std::array<unsigned, 2> M) : L(L), M(M), N(M[0] * M[1], 0), F(M[0] * M[1], 0.0) {}
+
+    void add(coordinate r, double x) {
+      unsigned ind = (unsigned)floor(M[1] * (r.y / L.y)) + M[1] *
+        (unsigned)floor(M[0] * (r.x / L.x));
+      N[ind]++;
+      F[ind] += x;
+    }
+
+};
+
 class ma : public hooks {
   // need:
   //  - interface for turning on and off specific measurements
@@ -55,8 +73,7 @@ class ma : public hooks {
     std::vector<uint64_t> cq;
 
     std::list<std::array<double, 2>> aG;
-    std::vector<double> cS;
-    std::vector<uint64_t> NS;
+    RealField fS;
     std::vector<uint64_t> sp;
 
     ClusterTree last_clusters;