Batch system can now plot ECGs

.gitignore

@@ -89,6 +89,7 @@ cuboid_tt2004_gpu/
 *.ctags
 tags
 scripts/output
+scripts/private_scripts
 scripts/animations/convert_png_to_mp4.sh
 scripts/animations/frames
 scripts/animations/drop_zone
@@ -184,6 +185,7 @@ sync.sh
 sync_local.sh
 sync_outputs.sh
 private_configs/
+private_parameter_sets/
 /src/domains_library/custom_mesh_info_data.h
 retired_code/
 *.gcno

build.sh

@@ -87,6 +87,7 @@ for i in "${BUILD_ARGS[@]}"; do
             COMPILE_FIBER_CONVERTER='y'
             COMPILE_POSTPROCESSOR='y'
             #COMPILE_EXPAND='y'
+            COMPILE_CLIP='y'
             ;;
         simulator)
             COMPILE_SIMULATOR='y'
@@ -271,6 +272,11 @@ if [ -n "$COMPILE_EXPAND" ]; then
     COMPILE_EXECUTABLE "MonoAlg3D_expand_scar" "src/main_expand_scar.c" "" "$STATIC_DEPS" "$DYNAMIC_DEPS" "$CUDA_LIBRARY_PATH $CUDA_MATH_LIBRARY_PATH $EXTRA_LIB_PATH $LIBRARY_OUTPUT_DIRECTORY"
 fi
 
+if [ -n "$COMPILE_CLIP" ]; then
+    COMPILE_EXECUTABLE "MonoAlg3D_clip_mesh" "src/main_clip_mesh.c" "" "$STATIC_DEPS" "$DYNAMIC_DEPS" "$CUDA_LIBRARY_PATH $CUDA_MATH_LIBRARY_PATH $EXTRA_LIB_PATH $LIBRARY_OUTPUT_DIRECTORY"
+fi
+
+
 FIND_CRITERION
 
 if [ -n "$CRITERION_FOUND" ]; then

scripts/conductivity_formula.py

@@ -0,0 +1,11 @@
+def calc_sigma (intra_conductivity, extra_conductivity):
+	return (intra_conductivity*extra_conductivity)/(intra_conductivity+extra_conductivity)
+
+# N-benchmark conductivities
+sigma_i_l = 0.17		# Intracellular longitudinal {S/m}
+sigma_e_l = 0.62		# Extracellular longitudinal {S/m}
+sigma_i_t = 0.019		# Intracellular transversal {S/m}
+sigma_e_t = 0.24		# Extracellular transversal {S/m}
+
+print("sigma_l = %g" % (calc_sigma(sigma_i_l, sigma_e_l)))
+print("sigma_t = %g" % (calc_sigma(sigma_i_t, sigma_e_t)))

src/config/ecg_config.h

@@ -14,16 +14,16 @@
 #define CALC_ECG(name) void name(struct time_info *time_info, struct config *config, struct grid *the_grid)
 typedef CALC_ECG(calc_ecg_fn);
 
-#define INIT_CALC_ECG(name) void name(struct config *config, struct ode_solver *the_ode_solver, struct grid *the_grid)
+#define INIT_CALC_ECG(name) void name(struct config *config, struct ode_solver *the_ode_solver, struct grid *the_grid, char *output_dir)
 typedef INIT_CALC_ECG(init_calc_ecg_fn);
 
 #define END_CALC_ECG(name) void name(struct config *config)
 typedef END_CALC_ECG(end_calc_ecg_fn);
 
-#define CALL_INIT_CALC_ECG(config, ode_solver, grid)                                                                                                           \
+#define CALL_INIT_CALC_ECG(config, ode_solver, grid, output_dir)                                                                                                           \
     do {                                                                                                                                                       \
         if((config) && (config)->init_function) {                                                                                                              \
-            ((init_calc_ecg_fn *)(config)->init_function)((config), (ode_solver), (grid));                                                                           \
+            ((init_calc_ecg_fn *)(config)->init_function)((config), (ode_solver), (grid), (output_dir));                                                                           \
         }                                                                                                                                                      \
     } while(0)
 

src/ecg_library/ecg.c

@@ -185,9 +185,11 @@ static void get_leads(struct config *config) {
 INIT_CALC_ECG(init_pseudo_bidomain_cpu) {
     config->persistent_data = CALLOC_ONE_TYPE(struct pseudo_bidomain_persistent_data);
 
-    char *filename = strdup("./ecg.txt");
-    GET_PARAMETER_STRING_VALUE_OR_USE_DEFAULT(filename, config, "filename");
-
+    // The filename for the ECG will be always the "OUTPUT_DIR/ecg.txt"
+    uint32_t nlen_output_dir = strlen(output_dir);
+    char *filename = (char*)malloc(sizeof(char)*nlen_output_dir+10);
+    sprintf(filename, "%s/ecg.txt", output_dir);
+
     char *dir = get_dir_from_path(filename);
     create_dir(dir);
     free(dir);
@@ -285,7 +287,7 @@ END_CALC_ECG(end_pseudo_bidomain_cpu) {
 
 INIT_CALC_ECG(init_pseudo_bidomain_gpu) {
 
-    init_pseudo_bidomain_cpu(config, NULL, the_grid);
+    init_pseudo_bidomain_cpu(config, NULL, the_grid, output_dir);
 
     if(!the_ode_solver->gpu) {
         log_warn("The current implementation of pseudo_bidomain_gpu only works when the odes are also being solved using the GPU! Falling back to CPU version\n");
@@ -465,13 +467,13 @@ INIT_CALC_ECG(init_pseudo_bidomain) {
     GET_PARAMETER_BOOLEAN_VALUE_OR_USE_DEFAULT(gpu, config, "use_gpu");
     if(gpu) {
 #ifdef COMPILE_CUDA
-        init_pseudo_bidomain_gpu(config, the_ode_solver, the_grid);
+        init_pseudo_bidomain_gpu(config, the_ode_solver, the_grid, output_dir);
 #else
         log_warn("Cuda runtime not found in this system. Falling back to CPU version!!\n");
-        init_pseudo_bidomain_cpu(config, NULL, the_grid);
+        init_pseudo_bidomain_cpu(config, NULL, the_grid, output_dir);
 #endif
     } else {
-        init_pseudo_bidomain_cpu(config, NULL, the_grid);
+        init_pseudo_bidomain_cpu(config, NULL, the_grid, output_dir);
     }
 }
 

src/main_clip_mesh.c

@@ -0,0 +1,193 @@
+// Author: Lucas Berg (@bergolho)
+// Script used to clip a section of a mesh and extract the extra data information
+// Version: 29/01/2024
+// Last change: 29/01/2024
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <getopt.h>
+
+#include "3dparty/stb_ds.h"
+#include "common_types/common_types.h"
+#include "utils/file_utils.h"
+#include "config/config_common.h"
+#include "config/domain_config.h"
+#include "domains_library/mesh_info_data.h"
+#include "extra_data_library/helper_functions.h"
+
+static const char *expansion_opt_string = "i:o:n:v:h?";
+static const struct option long_expansion_options[] = {{"input", required_argument, NULL, 'i'},
+                                                        {"output", required_argument, NULL, 'o'},
+                                                        {"n_rows", required_argument, NULL, 'n'},
+                                                        {"n_volumes", required_argument, NULL, 'v'}};
+struct expansion_options {
+    char *input;
+    char *output;
+    int n_rows;
+    char* n_volumes;
+};
+
+static void display_expansion_usage(char **argv) {
+
+    printf("Usage: %s [options]\n\n", argv[0]);
+    printf("Options:\n");
+    printf("--input  | -i [input]. File. Default NULL.\n");
+    printf("--output | -o [output]. Output file. Default NULL.\n");
+    printf("--n_rows | -n [n_rows]. Number of rows to expand.\n");
+    printf("--n_volumes | -v [n_volumes]. Number of volumes in the mesh\n");
+    printf("--help | -h. Shows this help and exit \n");
+    exit(EXIT_FAILURE);
+}
+
+
+static void parse_expansion_options(int argc, char **argv, struct expansion_options *user_args) {
+
+    int opt = 0;
+    int option_index;
+
+    opt = getopt_long_only(argc, argv, expansion_opt_string, long_expansion_options, &option_index);
+
+    while(opt != -1) {
+        switch(opt) {
+        case 'i':
+            user_args->input = strdup(optarg);
+            break;
+        case 'o':
+            user_args->output = strdup(optarg);
+            break;
+        case 'n':
+            user_args->n_rows = atoi(optarg);
+            break;
+        case 'v':
+            user_args->n_volumes = strdup(optarg);
+            break;
+        case 'h': /* fall-through is intentional */
+        case '?':
+            display_expansion_usage(argv);
+            break;
+        default:
+            /* You won't actually get here. */
+            break;
+        }
+
+        opt = getopt_long(argc, argv, expansion_opt_string, long_expansion_options, &option_index);
+    }
+
+    if(!user_args->input) {
+        display_expansion_usage(argv);
+    }
+}
+
+static void expand_scar(const char *input, const char *output, int n_rows, char* n_volumes) {
+
+    //set_no_stdout(true);
+
+    struct grid *grid = new_grid();
+    struct config *domain_config;
+
+    domain_config = alloc_and_init_config_data();
+    char *discretization = "300";
+    //char *discretization = "500";
+
+    shput_dup_value(domain_config->config_data, "start_discretization", strdup(discretization));
+    shput_dup_value(domain_config->config_data, "maximum_discretization", strdup(discretization));
+    shput_dup_value(domain_config->config_data, "mesh_file", strdup(input));
+
+    domain_config->main_function_name = strdup("initialize_grid_with_dti_mesh");
+    shput_dup_value(domain_config->config_data, "name", "Oxford DTI004 with Transmurality and Fiber orientation");
+
+    //shput(domain_config->config_data, "side_length_x", strdup(discretization));
+    //shput(domain_config->config_data, "side_length_y", strdup(discretization));
+    //shput(domain_config->config_data, "side_length_z", strdup(discretization));
+    shput(domain_config->config_data, "num_volumes", strdup(n_volumes));
+    //shput(domain_config->config_data, "num_extra_fields", "5");
+
+    init_config_functions(domain_config, "./shared_libs/libdefault_domains.so", "domain");
+
+    int success = ((set_spatial_domain_fn*)domain_config->main_function)(domain_config, grid);
+
+    if(success == 0) {
+        printf("Error loading mesh in %s. Exiting!\n", input);
+        exit(EXIT_FAILURE);
+    }
+
+    struct dti_mesh_info *extra_data = NULL;
+    bool ignore_cell;
+    real_cpu dx, dy, dz;
+    real_cpu min_x, max_x, min_y, max_y, min_z, max_z;
+    real *f, *s, *n;
+    real_cpu transmurality, base_apex_heterogeneity, apicobasal;
+    uint32_t transmurality_labels;
+
+    min_x=76803.6;
+    min_y=51176.3;
+    min_z=8740.9;
+    max_x=96803.6;
+    max_y=71176.3;
+    max_z=28740.9;
+
+    FILE *out = fopen(output, "w");
+    FOR_EACH_CELL(grid) {
+        if(cell->active) {
+
+            f = cell->sigma.fibers.f;
+            s = cell->sigma.fibers.s;
+            n = cell->sigma.fibers.n;
+
+            dx = cell->discretization.x / 2.0;
+            dy = cell->discretization.x / 2.0;
+            dz = cell->discretization.x / 2.0;
+
+            extra_data = (struct dti_mesh_info *)cell->mesh_extra_info;
+            transmurality = extra_data->transmurality;
+            base_apex_heterogeneity = extra_data->base_apex_heterogeneity;
+            apicobasal = extra_data->apicobasal;
+            transmurality_labels = extra_data->dti_transmurality_labels;
+
+            // Change the FAST_ENDO tag to ENDO
+            if (transmurality_labels == 3) {
+                transmurality_labels = 0;
+            }
+
+            ignore_cell = cell->center.x < min_x || cell->center.x > max_x || cell->center.y < min_y || cell->center.y > max_y || cell->center.z < min_z || cell->center.z > max_z;
+
+            if(!ignore_cell) {
+                fprintf(out, "%g,%g,%g,%g,%g,%g,%g,%g,%g,%u,%g,%g,%g,%g,%g,%g,%g,%g,%g\n", cell->center.x, cell->center.y, cell->center.z, dx, dy, dz, \
+                                                    transmurality, base_apex_heterogeneity, apicobasal, \
+                                                    transmurality_labels, \
+                                                    f[0], f[1], f[2], s[0], s[1], s[2], n[0], n[1], n[2]);
+            }
+
+        }
+    }
+    fclose(out);
+}
+
+int main(int argc, char **argv) {
+
+    struct expansion_options *options  = CALLOC_ONE_TYPE(struct expansion_options);
+
+    parse_expansion_options(argc, argv, options);
+
+    char *input = options->input;
+    char *output = options->output;
+
+    struct path_information input_info;
+
+    get_path_information(input, &input_info);
+
+    if(!input_info.exists) {
+        fprintf(stderr, "Invalid file (%s)! The input parameter should be an existing alg file!\n", input);
+        return EXIT_FAILURE;
+    }
+
+    if(!output) {
+        output = "expanded_mesh.alg";
+    }
+
+    expand_scar(input, output, options->n_rows, options->n_volumes);
+
+
+    return EXIT_SUCCESS;
+}

src/matrix_assembly_library/matrix_assembly.c

@@ -1220,79 +1220,3 @@ ASSEMBLY_MATRIX(anisotropic_sigma_assembly_matrix_with_fast_endocardium_layer) {
     free(s);
     free(n);
 }
-
-// Albert`s code
-ASSEMBLY_MATRIX(anisotropic_sigma_assembly_matrix_with_scaling) {
-
-    uint32_t num_active_cells = the_grid->num_active_cells;
-    struct cell_node **ac = the_grid->active_cells;
-
-    initialize_diagonal_elements(the_solver, the_grid);
-
-    real_cpu D[3][3];
-    int i;
-
-    real_cpu sigma_l = 0.0;
-    real_cpu sigma_t = 0.0;
-    real_cpu sigma_n = 0.0;
-
-    char *fiber_file = NULL;
-    GET_PARAMETER_STRING_VALUE_OR_USE_DEFAULT(fiber_file, config, "fibers_file");
-
-    struct fiber_coords_scale *fibers = NULL;
-    struct fiber_coords *fibers2 = NULL;
-
-    GET_PARAMETER_NUMERIC_VALUE_OR_REPORT_ERROR(real_cpu, sigma_l, config, "sigma_l");
-    GET_PARAMETER_NUMERIC_VALUE_OR_REPORT_ERROR(real_cpu, sigma_t, config, "sigma_t");
-    GET_PARAMETER_NUMERIC_VALUE_OR_REPORT_ERROR(real_cpu, sigma_n, config, "sigma_n");
-
-    real_cpu *f = NULL;
-    real_cpu *s = NULL;
-    real_cpu *n = NULL;
-    real_cpu *x = NULL;
-
-    if(fiber_file) {
-        log_info("Loading mesh fibers\n");
-        fibers = read_fibers_scale(fiber_file);
-        fibers2 = read_fibers(fiber_file, false);
-    }
-    else {
-        log_error_and_exit("Fibers file not provided");
-    }
-
-    OMP(parallel for private(D))
-    for(i = 0; i < num_active_cells; i++) {
-
-        int fiber_index = ac[i]->original_position_in_file;
-
-        if(fiber_index == -1) {
-            log_error_and_exit("fiber_index should not be -1, but it is for cell in index %d - %lf, %lf, %lf\n", i, ac[i]->center.x, ac[i]->center.y, ac[i]->center.z);
-        }
-
-        if(sigma_t == sigma_n) {
-            calc_tensor2(D, fibers[fiber_index].f, sigma_l*fibers[fiber_index].x[0], sigma_t*fibers[fiber_index].x[1]);
-        }
-        else {
-            calc_tensor(D, fibers[fiber_index].f, fibers[fiber_index].s, fibers[fiber_index].n, sigma_l*fibers[fiber_index].x[0], sigma_t*fibers[fiber_index].x[1], sigma_n*fibers[fiber_index].x[2]);
-        }
-        ac[i]->sigma.fibers = fibers2[fiber_index];
-
-        ac[i]->sigma.x = D[0][0];
-        ac[i]->sigma.y = D[1][1];
-        ac[i]->sigma.z = D[2][2];
-
-        ac[i]->sigma.xy = D[0][1];
-        ac[i]->sigma.xz = D[0][2];
-        ac[i]->sigma.yz = D[1][2];
-    }
-
-    OMP(parallel for)
-    for(i = 0; i < num_active_cells; i++) {
-        fill_discretization_matrix_elements_aniso(ac[i]);
-    }
-
-    free(f);
-    free(s);
-    free(n);
-    free(x);
-}

src/monodomain/monodomain_solver.c

@@ -601,7 +601,7 @@ int solve_monodomain(struct monodomain_solver *the_monodomain_solver, struct ode
 
     CALL_INIT_LINEAR_SYSTEM(linear_system_solver_config, the_grid, false || !domain_config);
     CALL_INIT_SAVE_MESH(save_mesh_config);
-    CALL_INIT_CALC_ECG(calc_ecg_config, the_ode_solver, the_grid);
+    CALL_INIT_CALC_ECG(calc_ecg_config, the_ode_solver, the_grid, out_dir_name);
 
     if(purkinje_linear_system_solver_config) {
         CALL_INIT_LINEAR_SYSTEM(purkinje_linear_system_solver_config, the_grid, true);