added energy dissipated

1 files changed, 53 insertions, 3 deletions

diff --git a/src/fracture.c b/src/fracture.c
index 79558c2..11e6efd 100644
--- a/src/fracture.c
+++ b/src/fracture.c
@@ -11,7 +11,7 @@ int main(int argc, char *argv[]) {
 	uint_t L;
 	double beta, inf, cutoff, crack_len;
 	bool crack_growth_crit, save_data, save_cluster_dist, use_voltage_boundaries, use_dual, save_network,
-			 save_crit_stress, save_stress_field, save_voltage_field, save_toughness, save_conductivity,
+			 save_crit_stress, save_stress_field, save_voltage_field, save_toughness, save_energy, save_conductivity,
 			 save_damage, save_damage_field, save_threshold;
 	bound_t boundary;
 	lattice_t lattice;
@@ -45,6 +45,7 @@ int main(int argc, char *argv[]) {
 	save_damage_field = false;
 	save_conductivity = false;
 	save_toughness = false;
+	save_energy = false;
 	save_threshold = false;
 
 
@@ -57,7 +58,7 @@ int main(int argc, char *argv[]) {
 
 	// get commandline options
 
-	while ((opt = getopt(argc, argv, "n:L:b:B:q:dVcoNsCrtDSvel:gT")) != -1) {
+	while ((opt = getopt(argc, argv, "n:L:b:B:q:dVcoNsCrtDSvel:gTE")) != -1) {
 		switch (opt) {
 			case 'n':
 				N = atoi(optarg);
@@ -154,6 +155,9 @@ int main(int argc, char *argv[]) {
 			case 't':
 				save_toughness = true;
 				break;
+			case 'E':
+				save_energy = true;
+				break;
 			case 'T':
 				save_threshold = true;
 				break;
@@ -265,6 +269,11 @@ int main(int argc, char *argv[]) {
 		toughness = (double *)malloc(N * sizeof(double));
 	}
 
+	double *energy;
+	if (save_energy) {
+		energy = (double *)malloc(N * sizeof(double));
+	}
+
 	double *thresholds;
 	if (save_threshold) {
 		thresholds = (double *)malloc(N * sizeof(double));
@@ -300,6 +309,13 @@ int main(int argc, char *argv[]) {
 		uint_t max_pos = 0;
 		double max_val = 0;
 
+		double cond0;
+		{
+			double *tmp_voltages = net_voltages(net, &c);
+			cond0 = net_conductivity(net, tmp_voltages);
+			free(tmp_voltages);
+		}
+
 		for (uint_t j = 0; j < data->num_broken; j++) {
 			double val = data->extern_field[j];
 
@@ -359,7 +375,13 @@ int main(int argc, char *argv[]) {
 
 		if (save_crit_stress) crit_stress[i] = data->extern_field[max_pos];
 
-		if (save_conductivity) conductivity[i] = data->conductivity[max_pos];
+		if (save_conductivity) {
+			if (max_pos > 0) {
+				conductivity[i] = data->conductivity[max_pos - 1];
+			} else {
+				conductivity[i] = cond0;
+			}
+		}
 
 		if (save_damage) damage[max_pos]++;
 
@@ -394,6 +416,24 @@ int main(int argc, char *argv[]) {
 			}
 		}
 
+		if (save_energy) {
+			double tmp_energy = 0;
+			if (max_pos > 0) {
+				double sigma1 = data->extern_field[0];
+				double cond1 = cond0;
+				for (uint_t j = 0; j < max_pos - 1; j++) {
+					double sigma2 = data->extern_field[j+1];
+					double cond2 = data->conductivity[j];
+					if (sigma2 > sigma1) {
+						tmp_energy += 0.5 * gsl_pow_2(sigma1) * (1 - cond2 / cond1) / cond1;
+						sigma1 = sigma2;
+						cond1 = cond2;
+					}
+				}
+			}
+			energy[i] = tmp_energy;
+		}
+
 		if (save_toughness) {
 			double tmp_toughness = 0;
 			if (max_pos > 0) {
@@ -531,6 +571,16 @@ int main(int argc, char *argv[]) {
 		free(toughness);
 	}
 
+	if (save_energy) {
+		char *tough_filename = (char *)malloc(filename_len * sizeof(char));
+		snprintf(tough_filename, filename_len, "enrg_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);
+		FILE *tough_file = fopen(tough_filename, "ab");
+		fwrite(energy, sizeof(double), N, tough_file);
+		fclose(tough_file);
+		free(tough_filename);
+		free(energy);
+	}
+
 	if (save_threshold) {
 		char *thres_filename = (char *)malloc(filename_len * sizeof(char));
 		snprintf(thres_filename, filename_len, "thrs_%c_%c_%c_%c_%d_%g_%g.dat", lattice_c, dual_c, boundc, boundc2, L, beta, crack_len);