Refactoring the Purkinje coupling code/Minor fixes in the conic domai…

…n function

ToDo

@@ -19,3 +19,12 @@
 The logger symbols are only exported to an executable if an static library linked to a shared library uses then. For now this is ok. But I think it will be a pain in future releases.
 The GPU linear system solver is not working for purkinje-only simulations
 When the minimum number of PMJs is not reached the solver will be in an infinite loop
+
+Lucas ToDo list:
+
+0) Revise all the examples inside the example_configs/ folder.
+
+1) Purkinje coupling: 
+    - Find a better way to write the PMJ delay into a file for the propagation block cases.
+
+

example_configs/cable_mesh_with_ToRORd_fkatp_endo_mid_epi.ini

@@ -1,6 +1,10 @@
 # ====================================================================
 # Author: Lucas Berg
 # Description: Simple simulation to test the ToRORd_fkatp model in a cable.
+# With the [extra_data] section we consider:
+#   45% of the cable ENDO
+#   25% of the cable MID
+#   30% of the cable EPI
 # When no [extra_data] information is provided all cells are 
 # considered to be control ENDO.
 # ====================================================================

networks/simple_his_bundle_pmj_location.vtk

@@ -0,0 +1,10 @@
+# vtk DataFile Version 4.2
+vtk output
+ASCII
+DATASET POLYDATA
+POINTS 2 float
+10000 15000 0
+10000 5000 0
+VERTICES 2 4
+1 0
+1 1

scripts/error_calculator/main.cpp

@@ -71,6 +71,10 @@ int main (int argc, char *argv[])
         cerr << "[-] ERROR! The number of cells are different" << endl;
         exit(EXIT_FAILURE);
     }
+    else
+    {
+        cout << "[+] Cell number are equal " << total_num_cells_1 << " " << total_num_cells_2 << endl;
+    }
 
     vtkSmartPointer<vtkFloatArray> array_1 = 
             vtkFloatArray::SafeDownCast(unstructured_grid_1->GetCellData()->GetArray(array_name.c_str()));

src/domains_library/domain.c

@@ -439,8 +439,11 @@ SET_SPATIAL_DOMAIN(initialize_grid_with_cuboid_and_sphere_fibrotic_mesh_with_con
     real_cpu phi = 0.0;
     GET_PARAMETER_NUMERIC_VALUE_OR_REPORT_ERROR(real_cpu, phi, config, "phi");
 
-    real_cpu plain_center = 0.0;
-    GET_PARAMETER_NUMERIC_VALUE_OR_REPORT_ERROR(real_cpu, plain_center, config, "plain_center");
+    real_cpu plain_center_x = 0.0;
+    GET_PARAMETER_NUMERIC_VALUE_OR_REPORT_ERROR(real_cpu, plain_center_x, config, "plain_center_x");
+
+    real_cpu plain_center_y = 0.0;
+    GET_PARAMETER_NUMERIC_VALUE_OR_REPORT_ERROR(real_cpu, plain_center_y, config, "plain_center_y");
 
     real_cpu sphere_radius = 0.0;
     GET_PARAMETER_NUMERIC_VALUE_OR_REPORT_ERROR(real_cpu, sphere_radius, config, "sphere_radius");
@@ -459,7 +462,7 @@ SET_SPATIAL_DOMAIN(initialize_grid_with_cuboid_and_sphere_fibrotic_mesh_with_con
 
     //set_square_mesh(config, the_grid);
     set_cuboid_domain_mesh(the_grid, start_dx, start_dy, start_dz, side_length_x, side_length_y, side_length_z);
-    set_cuboid_sphere_fibrosis_with_conic_path(the_grid, phi, plain_center, sphere_radius, border_zone_size, border_zone_radius, seed, conic_slope);
+    set_cuboid_sphere_fibrosis_with_conic_path(the_grid, phi, plain_center_x, plain_center_y, sphere_radius, border_zone_size, border_zone_radius, seed, conic_slope);
 
     return 1;
-}
+}

src/domains_library/domain_helpers.c

@@ -1010,8 +1010,8 @@ int calc_num_refs(real_cpu start_h, real_cpu desired_h) {
     return num_refs;
 }
 
-void set_cuboid_sphere_fibrosis_with_conic_path(struct grid *the_grid, real_cpu phi, real_cpu plain_center, real_cpu sphere_radius, real_cpu bz_size, real_cpu bz_radius,
-                               unsigned fib_seed, real_cpu cone_slope) {
+void set_cuboid_sphere_fibrosis_with_conic_path(struct grid *the_grid, real_cpu phi, real_cpu plain_center_x, real_cpu plain_center_y, \
+                                                real_cpu sphere_radius, real_cpu bz_size, real_cpu bz_radius, unsigned fib_seed, real_cpu cone_slope) {
 
     log_info("Making %.2lf %% of cells inactive\n", phi * 100.0f);
 
@@ -1028,23 +1028,22 @@ void set_cuboid_sphere_fibrosis_with_conic_path(struct grid *the_grid, real_cpu
 
     grid_cell = the_grid->first_cell;
     while(grid_cell != 0) {
-        //Calcula distância da célula para o centro da malha
-        real_cpu distance = pow(grid_cell->center.x - plain_center, 2.0) + pow(grid_cell->center.y - plain_center, 2.0);
-        real_cpu h_distance = abs(grid_cell->center.x - plain_center); 
+        // Calculate distance to the center of the mesh
+        real_cpu distance = pow(grid_cell->center.x - plain_center_x, 2.0) + pow(grid_cell->center.y - plain_center_y, 2.0);
+        real_cpu h_distance = abs(grid_cell->center.y - plain_center_y); 
 
         if(grid_cell->active) {
 
             INITIALIZE_FIBROTIC_INFO(grid_cell);
 
             if(distance <= bz_radius_2) {
-                //Dentro da border zone 
+                // Inside border zone 
                 if(distance <= sphere_radius_2) {
-                    //dentro da "esfera"
-                    if(h_distance < cone_slope*grid_cell->center.y*plain_center) //abs(cone_slope*grid_cell->center.y)
+                    // Inside the sphere
+                    if(h_distance < cone_slope*abs(the_grid->mesh_side_length.x - grid_cell->center.x)*plain_center_x) //abs(cone_slope*grid_cell->center.y)
                     {
+                        // Inside the cone
                         FIBROTIC(grid_cell) = true;
-
-                        //Dentro do cone
                     }
                     else{
                         grid_cell->active = false;
@@ -1069,8 +1068,8 @@ void set_cuboid_sphere_fibrosis_with_conic_path(struct grid *the_grid, real_cpu
                     grid_cell->active = false;
                 grid_cell->can_change = false;
             } else if(BORDER_ZONE(grid_cell)) {
-                real_cpu distance_from_center = sqrt((grid_cell->center.x - plain_center) * (grid_cell->center.x - plain_center) +
-                                                     (grid_cell->center.y - plain_center) * (grid_cell->center.y - plain_center));
+                real_cpu distance_from_center = sqrt((grid_cell->center.x - plain_center_x) * (grid_cell->center.x - plain_center_x) +
+                                                     (grid_cell->center.y - plain_center_y) * (grid_cell->center.y - plain_center_y));
                 distance_from_center = (distance_from_center - sphere_radius) / bz_size;
                 real_cpu phi_local = phi - phi * distance_from_center;
                 real_cpu p = (real_cpu)(rand()) / (RAND_MAX);

src/domains_library/domain_helpers.h

@@ -56,9 +56,9 @@ uint32_t set_custom_mesh_from_file(struct grid *the_grid, const char *mesh_file,
 
 void set_cube_sphere_fibrosis(struct grid *the_grid, real_cpu phi, real_cpu sphere_center[3], real_cpu sphere_radius, unsigned fib_seed);
 
-void set_cuboid_sphere_fibrosis_with_conic_path (struct grid *the_grid, real_cpu phi, real_cpu plain_center, real_cpu sphere_radius, real_cpu bz_size, real_cpu bz_radius,
-                               unsigned fib_seed, real_cpu cone_slope);
-
+void set_cuboid_sphere_fibrosis_with_conic_path(struct grid *the_grid, real_cpu phi, real_cpu plain_center_x, real_cpu plain_center_y, \
+                                                real_cpu sphere_radius, real_cpu bz_size, real_cpu bz_radius, unsigned fib_seed, real_cpu cone_slope);
+                                                
 int calc_num_refs(real_cpu start_h, real_cpu desired_h);
 
 #endif // MONOALG3D_DOMAIN_HELPERS_H

src/domains_library/mesh_info_data.h

@@ -10,6 +10,7 @@ struct default_fibrotic_mesh_info {
     int tissue_type;
 };
 
+// TODO: Move this struct and its macros to "custom_mesh_info_data.h"!
 struct dti_mesh_info {
 	enum dti_transmurality_labels {
 		DTI_FAST_ENDO,
@@ -20,6 +21,7 @@ struct dti_mesh_info {
 	float transmurality;
     float apicobasal;
 	float base_apex_heterogeneity;
+    float sf_IKs;
     int fast_endo;
 };
 
@@ -45,5 +47,6 @@ struct dti_mesh_info {
 #define DTI_MESH_BASE_APEX_HETEROGENEITY(grid_cell) (DTI_MESH_INFO(grid_cell))->base_apex_heterogeneity
 #define DTI_MESH_APICOBASAL(grid_cell) (DTI_MESH_INFO(grid_cell))->apicobasal
 #define DTI_MESH_FAST_ENDO(grid_cell) (DTI_MESH_INFO(grid_cell))->fast_endo
+#define DTI_MESH_SCALE_FACTOR_IKS(grid_cell) (DTI_MESH_INFO(grid_cell))->sf_IKs
 
 #endif /* __MESH_INFO_DATA_H */

src/extra_data_library/helper_functions.h

@@ -31,11 +31,16 @@ struct extra_data_for_tt3 {
     real GCaL_multiplicator;
     real INaCa_multiplicator;
     real Ikatp_multiplicator;
+    real GKr_multiplicator;
+    real GKs_multiplicator;
+    real alpha_multiplicator;
     real Vm_modifier;
     real *fibrosis;
     real *transmurality;
 };
 
+
+
 struct extra_data_for_torord {
     real INa_Multiplier; 
     real ICaL_Multiplier;
@@ -54,6 +59,8 @@ struct extra_data_for_torord {
     real IClb_Multiplier; 
     real Jrel_Multiplier; 
     real Jup_Multiplier;
+    real aCaMK_Multiplier;
+    real taurelp_Multiplier;
     real *initial_ss_endo;
     real *initial_ss_epi;
     real *initial_ss_mid;
@@ -120,11 +127,78 @@ struct extra_data_for_trovato {
     //real *purkinje_tags;
 };
 
+struct extra_data_for_torord_land_twave {
+    real INa_Multiplier; 
+    real INaL_Multiplier;
+    real INaCa_Multiplier;
+    real INaK_Multiplier;  
+    real INab_Multiplier;
+    real Ito_Multiplier;
+    real IKr_Multiplier; 
+    real IKs_Multiplier; 
+    real IK1_Multiplier; 
+    real IKb_Multiplier; 
+    real IKCa_Multiplier;
+    real ICaL_Multiplier;
+    real ICab_Multiplier;  
+    real IpCa_Multiplier;  
+    real ICaCl_Multiplier;
+    real IClb_Multiplier; 
+    real Jrel_Multiplier; 
+    real Jup_Multiplier;
+    real aCaMK_Multiplier;
+    real taurelp_Multiplier;
+    real *initial_ss_endo;
+    real *initial_ss_epi;
+    real *initial_ss_mid;
+    real *transmurality;
+    real *sf_IKs;
+    real *basetoapex;
+};
+
+struct extra_data_for_torord_general {
+    real INa_Multiplier; 
+    real INaL_Multiplier;
+    real INaCa_Multiplier;
+    real INaK_Multiplier;  
+    real INab_Multiplier;
+    real Ito_Multiplier;
+    real IKr_Multiplier; 
+    real IKs_Multiplier; 
+    real IK1_Multiplier; 
+    real IKb_Multiplier; 
+    real IKCa_Multiplier;
+    real ICaL_Multiplier;
+    real ICab_Multiplier;  
+    real IpCa_Multiplier;  
+    real ICaCl_Multiplier;
+    real IClb_Multiplier; 
+    real Jrel_Multiplier; 
+    real Jup_Multiplier;
+    real aCaMK_Multiplier;
+    real taurelp_Multiplier;
+    real *initial_ss_endo;
+    real *initial_ss_epi;
+    real *initial_ss_mid;
+    real *transmurality;
+    real *sf_IKs;
+    real *basetoapex;
+};
+
 struct extra_data_for_fibrosis * set_common_schemia_data(struct config *config, uint32_t num_cells);
 struct extra_data_for_tt3 * set_common_tt3_data (struct config *config, uint32_t num_cells);
 struct extra_data_for_torord * set_common_torord_data (struct config *config, uint32_t num_cells);
 struct extra_data_for_torord * set_common_torord_dyncl_data (struct config *config, uint32_t num_cells);
 struct extra_data_for_torord_land * set_common_torord_Land_data (struct config *config, uint32_t num_cells);
+struct extra_data_for_torord_land_twave * set_common_torord_Land_twave_data (struct config *config, uint32_t num_cells);
 struct extra_data_for_trovato * set_common_trovato_data (struct config *config, uint32_t num_cells);
 
+struct extra_data_for_torord_general * set_common_torord_general (struct config *config, uint32_t num_cells, uint32_t model_id);
+
+void configure_torord_extra_data_with_user_input (struct extra_data_for_torord_general *extra_data, struct config *config, uint32_t model_id);
+void configure_torord_extra_data_initial_conditions (struct extra_data_for_torord_general *extra_data, uint32_t model_id);
+void configure_torord_fkatp_initial_conditions (struct extra_data_for_torord_general *extra_data);
+void configure_torord_dynCl_initial_conditions (struct extra_data_for_torord_general *extra_data);
+void configure_torord_Land_initial_conditions (struct extra_data_for_torord_general *extra_data);
+
 #endif // MONOALG3D_C_EXTRA_DATA_HELPER_FUNCTIONS_H