diff --git a/.gitignore b/.gitignore index 02ba9797..9293c0ff 100644 --- a/.gitignore +++ b/.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 diff --git a/build.sh b/build.sh index aad24943..d9d6237c 100755 --- a/build.sh +++ b/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 diff --git a/parameter_sets/params_0.txt b/parameter_sets/params_0.txt new file mode 100644 index 00000000..c39a4710 --- /dev/null +++ b/parameter_sets/params_0.txt @@ -0,0 +1,8 @@ +2.0 +6.0 +138.3 +0.875 +0.875 +1.0 +1.7 +1.0 diff --git a/parameter_sets/params_1.txt b/parameter_sets/params_1.txt new file mode 100644 index 00000000..8efad49c --- /dev/null +++ b/parameter_sets/params_1.txt @@ -0,0 +1,8 @@ +2.0 +7.0 +138.3 +0.875 +0.875 +1.0 +1.7 +1.0 diff --git a/parameter_sets/params_2.txt b/parameter_sets/params_2.txt new file mode 100644 index 00000000..26a8f644 --- /dev/null +++ b/parameter_sets/params_2.txt @@ -0,0 +1,8 @@ +3.0 +6.0 +138.3 +0.875 +0.875 +1.0 +1.7 +1.0 diff --git a/parameter_sets/params_3.txt b/parameter_sets/params_3.txt new file mode 100644 index 00000000..a96c50b9 --- /dev/null +++ b/parameter_sets/params_3.txt @@ -0,0 +1,8 @@ +3.0 +7.0 +138.3 +0.875 +0.875 +1.0 +1.7 +1.0 diff --git a/scripts/conductivity_formula.py b/scripts/conductivity_formula.py new file mode 100644 index 00000000..66715db6 --- /dev/null +++ b/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))) diff --git a/src/config/ecg_config.h b/src/config/ecg_config.h index aa19a0b1..ba936312 100644 --- a/src/config/ecg_config.h +++ b/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) diff --git a/src/ecg_library/ecg.c b/src/ecg_library/ecg.c index f3c16751..095ad2ea 100644 --- a/src/ecg_library/ecg.c +++ b/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); } } diff --git a/src/main_clip_mesh.c b/src/main_clip_mesh.c new file mode 100644 index 00000000..77ce2143 --- /dev/null +++ b/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 +#include +#include +#include + +#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; +} \ No newline at end of file diff --git a/src/matrix_assembly_library/matrix_assembly.c b/src/matrix_assembly_library/matrix_assembly.c index 0bfebed4..e697a79f 100644 --- a/src/matrix_assembly_library/matrix_assembly.c +++ b/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); -} \ No newline at end of file diff --git a/src/monodomain/monodomain_solver.c b/src/monodomain/monodomain_solver.c index 06d9cf97..6f2a3e37 100644 --- a/src/monodomain/monodomain_solver.c +++ b/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);