fixed the crack stress, adding new measurements for mean ordinary stress

2 files changed, 59 insertions, 29 deletions

diff --git a/src/measurements.cpp b/src/measurements.cpp
index e7aa528..446d29c 100644
--- a/src/measurements.cpp
+++ b/src/measurements.cpp
@@ -171,7 +171,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(1 + 2 * (unsigned)ceil(sqrt(n)), 2)), last_clusters(2 * n) {
+      cS(pow(4 * (unsigned)ceil(sqrt(n)), 2) / 2), NS(pow(4 * (unsigned)ceil(sqrt(n)), 2) / 2), cT(pow(4 * (unsigned)ceil(sqrt(n)), 2) / 2), NT(pow(4 * (unsigned)ceil(sqrt(n)), 2) / 2),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) + "_";
@@ -190,7 +190,6 @@ ma::ma(unsigned n, double a, double beta, double weight)
   Na = 0;
   NA = 0;
   Nq = 0;
-  NS = 0;
 }
 
 ma::ma(unsigned Lx, unsigned Ly, double beta, double weight)
@@ -200,7 +199,7 @@ 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 / 2 + 1) * (Ly / 2 + 1)), last_clusters(Lx * Ly / 2) {
+      cq((Lx / 2 + 1) * (Ly / 2 + 1)), cS(Lx * Ly / 2), NS(Lx * Ly / 2), cT(Lx * Ly / 2), NT(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) + "_";
@@ -219,7 +218,6 @@ ma::ma(unsigned Lx, unsigned Ly, double beta, double weight)
   Na = 0;
   NA = 0;
   Nq = 0;
-  NS = 0;
 }
 
 ma::~ma() {
@@ -244,7 +242,10 @@ 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", NS, cS, model_string);
+  update_field_file("cS", 0, cS, model_string);
+  update_field_file("NS", 0, NS, model_string);
+  update_field_file("cT", 0, cT, model_string);
+  update_field_file("NT", 0, NT, model_string);
 }
 
 void ma::pre_fracture(const network&) {
@@ -485,14 +486,6 @@ void ma::post_fracture(network& n) {
       right_end = end_0e;
     }
 
-    std::cout << crack.first[0] << " " << crack.first[1] << "\n";
-    std::cout << left_end << " " << right_end << "\n";
-    for (unsigned e : biggest_fragment) {
-      std::cout << e << " ";
-    }
-    std::cout << "\n";
-    std::cout << Δr.x << " " << Δr.y << "\n";
-
     for (unsigned i = 0; i < n.G.dual_edges.size(); i++) {
       if (!last_backbone[i]) {
         double Δx_tmpl, Δy_tmpl, Δx_tmpr, Δy_tmpr, Δxr, Δxl, Δyr, Δyl;
@@ -504,10 +497,10 @@ void ma::post_fracture(network& n) {
           Δy_tmpl = fabs(n.G.dual_edges[left_end].r.y - n.G.dual_edges[i].r.y);
           Δyl = Δy_tmpl > n.G.L.y / 2 ? n.G.L.y - Δy_tmpl : Δy_tmpl;
 
-          Δx_tmpr = n.G.dual_edges[i].r.x - n.G.dual_edges[left_end].r.x;
+          Δx_tmpr = n.G.dual_edges[i].r.x - n.G.dual_edges[right_end].r.x;
           Δxr = Δx_tmpr < 0 ? Δx_tmpr + n.G.L.x : Δx_tmpr;
 
-          Δy_tmpr = fabs(n.G.dual_edges[left_end].r.y - n.G.dual_edges[i].r.y);
+          Δy_tmpr = fabs(n.G.dual_edges[right_end].r.y - n.G.dual_edges[i].r.y);
           Δyr = Δy_tmpr > n.G.L.y / 2 ? n.G.L.y - Δy_tmpr : Δy_tmpr;
         } else {
           Δx_tmpl = n.G.dual_edges[left_end].r.y - n.G.dual_edges[i].r.y;
@@ -516,31 +509,66 @@ void ma::post_fracture(network& n) {
           Δy_tmpl = fabs(n.G.dual_edges[left_end].r.x - n.G.dual_edges[i].r.x);
           Δyl = Δy_tmpl > n.G.L.x / 2 ? n.G.L.x - Δy_tmpl : Δy_tmpl;
 
-          Δx_tmpr = n.G.dual_edges[i].r.y - n.G.dual_edges[left_end].r.y;
+          Δx_tmpr = n.G.dual_edges[i].r.y - n.G.dual_edges[right_end].r.y;
           Δxr = Δx_tmpr < 0 ? Δx_tmpr + n.G.L.y : Δx_tmpr;
 
-          Δy_tmpr = fabs(n.G.dual_edges[left_end].r.x - n.G.dual_edges[i].r.x);
+          Δy_tmpr = fabs(n.G.dual_edges[right_end].r.x - n.G.dual_edges[i].r.x);
           Δyr = Δy_tmpr > n.G.L.x / 2 ? n.G.L.x - Δy_tmpr : Δy_tmpr;
         }
 
-        if (Δxl < Δxr) {
-          cS[(unsigned)floor((1 + 2 * (Mx / 2)) * (Δyl / n.G.L.y)) +
-             (Mx / 2 + 1) * (unsigned)floor((1 + 2 * (My / 2)) * (Δxl / n.G.L.x))] +=
-              last_currents.currents[i] / last_currents.conductivity[crack.first[0] % 2];
-        } else {
-          cS[(unsigned)floor((1 + 2 * (Mx / 2)) * (Δyr / n.G.L.y)) +
-             (Mx / 2 + 1) * (unsigned)floor((1 + 2 * (My / 2)) * (Δxr / n.G.L.x))] +=
-              last_currents.currents[i] / last_currents.conductivity[crack.first[0] % 2];
-        }
+        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];
       }
     }
 
-    NS++;
-
   } else {
     sp[0]++;
   }
 
+  for (unsigned i = 0; i < n.G.dual_edges.size(); i++) {
+    if (!last_backbone[i]) {
+      for (unsigned j = 0; j < n.G.dual_edges.size(); j++) {
+        if (!last_backbone[j]) {
+          double Δx_tmpl, Δy_tmpl, Δx_tmpr, Δy_tmpr, Δxr, Δxl, Δyr, Δyl;
+
+          if (crack.first[0] % 2 == 1) {
+            Δx_tmpl = fabs(n.G.dual_edges[j].r.x - n.G.dual_edges[i].r.x);
+            Δxl = Δx_tmpl > n.G.L.x / 2 ? n.G.L.x - Δx_tmpl : Δx_tmpl;
+
+            Δy_tmpl = fabs(n.G.dual_edges[j].r.y - n.G.dual_edges[i].r.y);
+            Δyl = Δy_tmpl > n.G.L.y / 2 ? n.G.L.y - Δy_tmpl : Δy_tmpl;
+          } else {
+            Δx_tmpl = fabs(n.G.dual_edges[j].r.y - n.G.dual_edges[i].r.y);
+            Δxl = Δx_tmpl > n.G.L.y / 2 ? n.G.L.y - Δx_tmpl : Δx_tmpl;
+
+            Δy_tmpl = fabs(n.G.dual_edges[j].r.x - n.G.dual_edges[i].r.x);
+            Δyl = Δy_tmpl > n.G.L.x / 2 ? n.G.L.x - Δy_tmpl : Δy_tmpl;
+
+          }
+
+          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];
+        }
+      }
+    }
+  }
+
   aG.push_back(last_currents.conductivity);
 }
 
diff --git a/src/measurements.hpp b/src/measurements.hpp
index c0458b5..0786cf8 100644
--- a/src/measurements.hpp
+++ b/src/measurements.hpp
@@ -29,7 +29,6 @@ class ma : public hooks {
     uint64_t Na;
     uint64_t NA;
     uint64_t Nq;
-    uint64_t NS;
     
     // measurement storage
     std::vector<uint64_t> sc;
@@ -57,6 +56,9 @@ class ma : public hooks {
 
     std::list<std::array<double, 2>> aG;
     std::vector<double> cS;
+    std::vector<uint64_t> NS;
+    std::vector<double> cT;
+    std::vector<uint64_t> NT;
     std::vector<uint64_t> sp;
 
     ClusterTree last_clusters;