New domains function for source-sink mismatch with fibrosis

example_configs/plain_mesh_source_sink_mismatch_with_fibrosis.ini

@@ -0,0 +1,158 @@
+[main]
+num_threads=6
+dt_pde=0.02
+simulation_time=5700.0
+abort_on_no_activity=false
+
+[update_monodomain]
+main_function=update_monodomain_default
+
+[save_result]
+print_rate = 1000
+output_dir = ./outputs/plain_mesh_source_sink_mismatch_with_fibrosis
+add_timestamp = false
+binary = true
+save_ode_state_variables = false
+main_function = save_as_ensight
+remove_older_simulation=true
+
+[assembly_matrix]
+init_function=set_initial_conditions_fvm
+sigma_x=0.00001334
+sigma_y=0.00001334
+sigma_z=0.00001334
+library_file=shared_libs/libdefault_matrix_assembly.so
+main_function=homogeneous_sigma_assembly_matrix
+
+[linear_system_solver]
+tolerance=1e-16
+use_preconditioner=no
+max_iterations=200
+library_file=shared_libs/libdefault_linear_system_solver.so
+use_gpu=yes
+main_function=conjugate_gradient
+init_function=init_conjugate_gradient
+end_function=end_conjugate_gradient
+
+[domain]
+name=Plain Mesh With Square Fibrotic Region
+num_layers=1
+start_dx=200.0
+start_dy=200.0
+start_dz=200.0
+side_length=50000.0
+phi=1.0
+seed=123456789
+region_min_x=8000.0
+region_max_x=42000.0
+region_min_y=8000.0
+region_max_y=42000.0
+region_min_z=0.0
+region_max_z=1000.0
+source_sink_min_x = 24000.0
+source_sink_max_x = 26000.0
+main_function=initialize_grid_with_square_mesh_and_source_sink_fibrotic_region
+
+[ode_solver]
+dt=0.02
+use_gpu=yes
+library_file=shared_libs/libten_tusscher_2006.so
+
+[stim_s1]
+start = 0.0
+duration = 2.0
+current = -50.0
+y_limit = 200.0
+main_function=stim_if_y_less_than
+
+[stim_s2]
+start = 500.0
+duration = 2.0
+current = -50.0
+y_limit = 200.0
+main_function=stim_if_y_less_than
+
+[stim_s3]
+start = 1000.0
+duration = 2.0
+current = -50.0
+y_limit = 200.0
+main_function=stim_if_y_less_than
+
+[stim_s4]
+start = 1450.0
+duration = 2.0
+current = -50.0
+y_limit = 200.0
+main_function=stim_if_y_less_than
+
+[stim_s5]
+start = 1900.0
+duration = 2.0
+current = -50.0
+y_limit = 200.0
+main_function=stim_if_y_less_than
+
+[stim_s6]
+start = 2325.0
+duration = 2.0
+current = -50.0
+y_limit = 200.0
+main_function=stim_if_y_less_than
+
+[stim_s7]
+start = 2750.0
+duration = 2.0
+current = -50.0
+y_limit = 200.0
+main_function=stim_if_y_less_than
+
+[stim_s8]
+start = 3150.0
+duration = 2.0
+current = -50.0
+y_limit = 200.0
+main_function=stim_if_y_less_than
+
+[stim_s9]
+start = 3550.0
+duration = 2.0
+current = -50.0
+y_limit = 200.0
+main_function=stim_if_y_less_than
+
+[stim_s10]
+start = 3925.0
+duration = 2.0
+current = -50.0
+y_limit = 200.0
+main_function=stim_if_y_less_than
+
+[stim_s11]
+start = 4300.0
+duration = 2.0
+current = -50.0
+y_limit = 200.0
+main_function=stim_if_y_less_than
+
+[stim_s12]
+start = 4650.0
+duration = 2.0
+current = -50.0
+y_limit = 200.0
+main_function=stim_if_y_less_than
+
+[stim_s13]
+start = 5000.0
+duration = 2.0
+current = -50.0
+y_limit = 200.0
+main_function=stim_if_y_less_than
+
+[stim_s14]
+start = 5350.0
+duration = 2.0
+current = -50.0
+y_limit = 200.0
+main_function=stim_if_y_less_than
+

src/domains_library/domain.c

@@ -361,3 +361,45 @@ SET_SPATIAL_DOMAIN(initialize_grid_with_custom_mesh) {
     return ret;
 
 }
+
+SET_SPATIAL_DOMAIN(initialize_grid_with_square_mesh_and_source_sink_fibrotic_region) {
+
+    real_cpu phi = 0.0;
+    GET_PARAMETER_NUMERIC_VALUE_OR_REPORT_ERROR(real_cpu, phi, config, "phi");
+
+    unsigned seed = 0;
+    GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(unsigned, seed, config, "seed");
+
+    real_cpu min_x = 0;
+    GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real_cpu, min_x, config, "region_min_x");
+
+    real_cpu max_x = 0;
+    GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real_cpu, max_x, config, "region_max_x");
+
+    real_cpu min_y = 0;
+    GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real_cpu, min_y, config, "region_min_y");
+
+    real_cpu max_y = 0;
+    GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real_cpu, max_y, config, "region_max_y");
+
+    real_cpu min_z = 0;
+    GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real_cpu, min_z, config, "region_min_z");
+
+    real_cpu max_z = 0;
+    GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real_cpu, max_z, config, "region_max_z");
+
+    real_cpu source_sink_min_x = 0.0;
+    GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real_cpu, source_sink_min_x, config, "source_sink_min_x");
+
+    real_cpu source_sink_max_x = 0.0;
+    GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real_cpu, source_sink_max_x, config, "source_sink_max_x");
+
+    real_cpu side_length = 0.0;
+    GET_PARAMETER_NUMERIC_VALUE_OR_REPORT_ERROR(real_cpu, side_length, config, "side_length");
+
+    set_square_mesh(config, the_grid);
+
+    set_plain_fibrosis_source_sink_region(the_grid, phi, seed, min_x, max_x, min_y, max_y, min_z, max_z, source_sink_min_x, source_sink_max_x, side_length);
+
+    return 1;
+}

src/domains_library/domain_helpers.c

@@ -1009,3 +1009,45 @@ int calc_num_refs(real_cpu start_h, real_cpu desired_h) {
 
     return num_refs;
 }
+
+void set_plain_fibrosis_source_sink_region (struct grid *the_grid, real_cpu phi, unsigned fib_seed, const double min_x, const double max_x, const double min_y,
+                                      const double max_y, const double min_z, const double max_z,
+                                      real_cpu source_sink_min_x, real_cpu source_sink_max_x, real_cpu side_length) {
+    log_info("Making %.2lf %% of cells inside the region inactive\n", phi * 100.0);
+
+    struct cell_node *grid_cell;
+
+    if(fib_seed == 0)
+        fib_seed = (unsigned)time(NULL) + getpid();
+
+    srand(fib_seed);
+
+    log_info("Using %u as seed\n", fib_seed);
+
+    real_cpu a1 = (2.0*side_length) / (side_length - 2*source_sink_min_x);
+    real_cpu b1 = -source_sink_min_x*a1;
+    real_cpu a2 = (2.0*side_length) / (side_length - 2*source_sink_max_x);
+    real_cpu b2 = -source_sink_max_x*a2;
+
+    grid_cell = the_grid->first_cell;
+    while(grid_cell != 0) {
+        real center_x = grid_cell->center.x;
+        real center_y = grid_cell->center.y;
+        real center_z = grid_cell->center.z;
+
+        if(center_x >= min_x && center_x <= max_x && center_y >= min_y && center_y <= max_y && center_z >= min_z && center_z <= max_z
+            && (center_y > a1*center_x + b1 || center_y > a2*center_x + b2)) {
+            if(grid_cell->active) {
+                real_cpu p = (real_cpu)(rand()) / (RAND_MAX);
+                if(p < phi) {
+                    grid_cell->active = false;
+                }
+
+                INITIALIZE_FIBROTIC_INFO(grid_cell);
+                FIBROTIC(grid_cell) = true;
+            }
+        }
+
+        grid_cell = grid_cell->next;
+    }
+}

src/domains_library/domain_helpers.h

@@ -46,6 +46,10 @@ void set_plain_fibrosis_using_file(struct grid *the_grid, const char filename[])
 void set_plain_fibrosis_inside_region(struct grid *the_grid, real_cpu phi, unsigned fib_seed, double min_x, double max_x, double min_y, double max_y,
                                       double min_z, double max_z);
 
+void set_plain_fibrosis_source_sink_region (struct grid *the_grid, real_cpu phi, unsigned fib_seed, const double min_x, const double max_x, const double min_y,
+                                      const double max_y, const double min_z, const double max_z,
+                                      real_cpu source_sink_min_x, real_cpu source_sink_max_x, real_cpu side_length);
+
 uint32_t set_custom_mesh_from_file(struct grid *the_grid, const char *mesh_file, uint32_t num_volumes, double start_h, uint8_t num_extra_fields,
                                    set_custom_data_for_mesh_fn set_custom_data_for_mesh);