Preparing repository for pull request to 'rsachetto-master' ...

Author: bergolho

ToDo

@@ -18,4 +18,4 @@
 #KNOW ISSUES:
 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  
+When the minimum number of PMJs is not reached the solver will be in an infinite loop

example_configs/cable_mesh_with_ToRORd_fkatp_endo_mid_epi.ini

@@ -1,10 +1,6 @@
 # ====================================================================
 # 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.
 # ====================================================================
@@ -80,4 +76,4 @@ x_limit = 500.0
 main_function=stim_if_x_less_than
 
 [extra_data]
-main_function=set_extra_data_mixed_torord_fkatp_epi_mid_endo
+main_function=set_extra_data_mixed_torord_fkatp_epi_mid_endo

scripts/evaluateBenchmarkCV/main.cpp

@@ -1,4 +1,7 @@
-// Author: Lucas Berg
+// ------------------------------------------------------------------------------------
+// Author: Lucas Berg, Julia Camps and Jenny Wang
+// Script to evaluate the conductivities calibrated using the 'tuneCV' on a 3D wedge.
+// ------------------------------------------------------------------------------------
 
 #include <iostream>
 #include <string>

scripts/tuneCV/main.cpp

@@ -28,7 +28,7 @@ using namespace std;
 const char MONOALG_PATH[500] = "/home/berg/Github/MonoAlg3D_C";
 // ----------------------------------------------------------
 
-const double TOLERANCE = 1.0e-03;
+const double TOLERANCE = 5.0e-03;
 
 double calculate_conduction_velocity_from_simulation ()
 {
@@ -135,6 +135,16 @@ void write_configuration_file (const double sigma) {
     fprintf(file,"library_file=%s/shared_libs/libdefault_matrix_assembly.so\n", MONOALG_PATH);
     fprintf(file,"main_function=homogeneous_sigma_assembly_matrix\n");
     fprintf(file,"\n");
+
+    //fprintf(file,"[assembly_matrix]\n");
+    //fprintf(file,"init_function=set_initial_conditions_fvm\n");
+    //fprintf(file,"sigma_l=%g\n",sigma);
+    //fprintf(file,"sigma_t=%g\n",sigma);
+    //fprintf(file,"sigma_n=%g\n",sigma);
+    //fprintf(file,"f=[1,0,0]\n");
+    //fprintf(file,"library_file=%s/shared_libs/libdefault_matrix_assembly.so\n", MONOALG_PATH);
+    //fprintf(file,"main_function=anisotropic_sigma_assembly_matrix\n");
+    //fprintf(file,"\n");
 
     fprintf(file,"[linear_system_solver]\n");
     fprintf(file,"tolerance=1e-16\n");

scripts/tuneCVbenchmark/CMakeLists.txt

@@ -1,12 +1,12 @@
 cmake_minimum_required(VERSION 2.8)
 
-PROJECT(tuneCV)
+PROJECT(tuneCVbenchmark)
 
 find_package(VTK REQUIRED)
 include(${VTK_USE_FILE})
 
 set( CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/bin )
 
-add_executable(tuneCV main.cpp )
+add_executable(tuneCVbenchmark main.cpp )
 
-target_link_libraries(tuneCV ${VTK_LIBRARIES})
+target_link_libraries(tuneCVbenchmark ${VTK_LIBRARIES})

scripts/tuneCVbenchmark/main.cpp

@@ -1,4 +1,7 @@
-// Author: Lucas Berg
+// -------------------------------------------------------------------------------------
+// Authors: Lucas Berg, Julia Camps and Jenny Wang
+// Script to calibrate the monodomain conductivities using a cable simulation.
+// -------------------------------------------------------------------------------------
 
 #include <iostream>
 #include <string>
@@ -29,8 +32,13 @@ using namespace std;
 
 const double TOLERANCE = 1.0e-02; // 1 cm/s
 
-double calculate_conduction_velocity_from_cable_simulation ()
-{
+// Change your MonoAlg3D path here:
+// ----------------------------------------------------------
+const char MONOALG_PATH[500] = "/home/berg/Github/MonoAlg3D_C";
+// ----------------------------------------------------------
+
+double calculate_conduction_velocity_from_cable_simulation () {
+    
     string filename = "outputs/cable/tissue_activation_time_map_pulse_it_0.vtu";
 
     // Read all the data from the file
@@ -94,33 +102,35 @@ double calculate_conduction_velocity_from_cable_simulation ()
 }
 
 // TODO: Maybe pass a pre-configured config file as an input parameter with the cellular model setup that the user will use
-void write_configuration_file (const double sigma)
-{
-    FILE *file = fopen("/home/jenny/MonoAlg3D_C/scripts/tuneCVbenchmark/configs/cable.ini","w+");
+void write_configuration_file (const double sigma) {
+    
+    char filename[500];
+    sprintf(filename,"%s/scripts/tuneCVbenchmark/configs/cable.ini",MONOALG_PATH);
+    FILE *file = fopen(filename,"w+");
 
     fprintf(file,"[main]\n");
     fprintf(file,"num_threads=6\n");
     fprintf(file,"dt_pde=0.01\n");
-    fprintf(file,"simulation_time=100.0\n");
+    fprintf(file,"simulation_time=150.0\n");
     fprintf(file,"abort_on_no_activity=false\n");
     fprintf(file,"use_adaptivity=false\n");
     fprintf(file,"quiet=true\n");
     fprintf(file,"\n");
 
     fprintf(file,"[update_monodomain]\n");
     fprintf(file,"main_function=update_monodomain_default\n");
-    fprintf(file,"library_file=/home/jenny/MonoAlg3D_C/shared_libs/libdefault_update_monodomain.so\n");
+    fprintf(file,"library_file=%s/shared_libs/libdefault_update_monodomain.so\n",MONOALG_PATH);
     fprintf(file,"\n");
 
     // For saving the LATs in a format that can be read for calculating the CVs
     fprintf(file,"[save_result]\n");
     fprintf(file,"print_rate=1\n");
-    fprintf(file,"output_dir=/home/jenny/MonoAlg3D_C/scripts/tuneCVbenchmark/outputs/cable\n");
+    fprintf(file,"output_dir=%s/scripts/tuneCVbenchmark/outputs/cable\n",MONOALG_PATH);
     fprintf(file,"save_pvd=true\n");
     fprintf(file,"file_prefix=V\n");
     fprintf(file,"save_activation_time=true\n");
     fprintf(file,"save_apd=false\n");
-    fprintf(file,"library_file=/home/jenny/MonoAlg3D_C/shared_libs/libdefault_save_mesh_purkinje.so\n");
+    fprintf(file,"library_file=%s/shared_libs/libdefault_save_mesh_purkinje.so\n",MONOALG_PATH);
     fprintf(file,"main_function=save_tissue_with_activation_times\n");
     fprintf(file,"init_function=init_save_tissue_with_activation_times\n");
     fprintf(file,"end_function=end_save_tissue_with_activation_times\n");
@@ -132,7 +142,7 @@ void write_configuration_file (const double sigma)
     fprintf(file,"sigma_x=%g\n",sigma);
     fprintf(file,"sigma_y=%g\n",sigma);
     fprintf(file,"sigma_z=%g\n",sigma);
-    fprintf(file,"library_file=/home/jenny/MonoAlg3D_C/shared_libs/libdefault_matrix_assembly.so\n");
+    fprintf(file,"library_file=%s/shared_libs/libdefault_matrix_assembly.so\n",MONOALG_PATH);
     fprintf(file,"main_function=homogeneous_sigma_assembly_matrix\n");
     fprintf(file,"\n");
 
@@ -141,7 +151,7 @@ void write_configuration_file (const double sigma)
     fprintf(file,"use_preconditioner=no\n");
     fprintf(file,"use_gpu=yes\n");
     fprintf(file,"max_iterations=200\n");
-    fprintf(file,"library_file=/home/jenny/MonoAlg3D_C/shared_libs/libdefault_linear_system_solver.so\n");
+    fprintf(file,"library_file=%s/shared_libs/libdefault_linear_system_solver.so\n",MONOALG_PATH);
     fprintf(file,"init_function=init_conjugate_gradient\n");
     fprintf(file,"end_function=end_conjugate_gradient\n");
     fprintf(file,"main_function=conjugate_gradient\n");
@@ -153,15 +163,15 @@ void write_configuration_file (const double sigma)
     fprintf(file,"start_dy=500.0\n");
     fprintf(file,"start_dz=500.0\n");
     fprintf(file,"cable_length=20000.0\n");
-    fprintf(file,"library_file=/home/jenny/MonoAlg3D_C/shared_libs/libdefault_domains.so\n");
+    fprintf(file,"library_file=%s/shared_libs/libdefault_domains.so\n",MONOALG_PATH);
     fprintf(file,"main_function=initialize_grid_with_cable_mesh\n");
     fprintf(file,"\n");
 
     fprintf(file,"[ode_solver]\n");
     fprintf(file,"dt=0.01\n");
     fprintf(file,"use_gpu=yes\n");
     fprintf(file,"gpu_id=0\n");
-    fprintf(file,"library_file=/home/jenny/MonoAlg3D_C/shared_libs/libten_tusscher_tt3_mixed_endo_mid_epi.so\n");
+    fprintf(file,"library_file=%s/shared_libs/libten_tusscher_tt3_mixed_endo_mid_epi.so\n",MONOALG_PATH);
     fprintf(file,"\n");
 
     fprintf(file,"[stim_benchmark]\n");
@@ -175,16 +185,14 @@ void write_configuration_file (const double sigma)
     fprintf(file, "min_z = 0.0\n");
     fprintf(file, "max_z = 3000.0\n");
     fprintf(file,"main_function=stim_x_y_z_limits\n");
-    fprintf(file,"library_file=/home/jenny/MonoAlg3D_C/shared_libs/libdefault_stimuli.so\n");
+    fprintf(file,"library_file=%s/shared_libs/libdefault_stimuli.so\n",MONOALG_PATH);
     fprintf(file,"\n");
 
     fclose(file);
 }
 
-int main (int argc, char *argv[])
-{
-    if (argc-1 != 1)
-    {
+int main (int argc, char *argv[]) {
+    if (argc-1 != 1) {
         cerr << "=============================================================================" << endl;
         cerr << "Usage:> " << argv[0] << " <target_CV>" << endl;
         cerr << "=============================================================================" << endl;
@@ -204,20 +212,21 @@ int main (int argc, char *argv[])
     double target_cv = atof(argv[1]);
     double sigma = 0.0002;
 
-    do
-    {
+    do {
         write_configuration_file(sigma);
 
         // Run the simulation
-        system("/home/jenny/MonoAlg3D_C/bin/MonoAlg3D -c /home/jenny/MonoAlg3D_C/scripts/tuneCVbenchmark/configs/cable.ini");
+        char command[500];
+        sprintf(command,"%s/bin/MonoAlg3D -c %s/scripts/tuneCVbenchmark/configs/cable.ini",MONOALG_PATH,MONOALG_PATH);
+        system(command);
 
         cv = calculate_conduction_velocity_from_cable_simulation();
         factor = pow(target_cv/cv,2);
         sigma = sigma*factor;
 
         printf("\n|| Target CV = %g m/s || Computed CV = %g m/s || Factor = %g || Adjusted sigma = %g mS/um ||\n\n",target_cv,cv,factor,sigma);
 
-    }while ( fabs(cv-target_cv) > TOLERANCE );
+    } while ( fabs(cv-target_cv) > TOLERANCE );
 
     printf("\n[+] Target conductivity = %g mS/um\n",sigma);
 

src/ensight_utils/ensight_grid.c

@@ -108,7 +108,7 @@ void save_case_file(char *filename, uint64_t num_files, real_cpu dt, int print_r
     fclose(case_file);
 }
 
-void save_en6_result_file(char *filename, struct grid *the_grid, bool binary) {
+void save_en6_result_file(char *filename, struct grid *the_grid, bool binary, bool save_purkinje) {
 
     FILE *result_file;
 
@@ -142,7 +142,7 @@ void save_en6_result_file(char *filename, struct grid *the_grid, bool binary) {
         part_number++;
     }
 
-    if(the_grid->purkinje) {
+    if(the_grid->purkinje && save_purkinje) {
         write_string("part", result_file, binary);
         new_line(result_file, binary);
 
@@ -366,26 +366,35 @@ static inline void set_point_data(struct point_3d center, struct point_3d half_f
 struct ensight_grid * new_ensight_grid_from_alg_grid(struct grid *grid, bool clip_with_plain,
                                                                      float *plain_coordinates, bool clip_with_bounds,
                                                                      float *bounds, bool read_fibers_f,
-                                                                     bool save_fibrotic) {
+                                                                     bool save_fibrotic, bool save_purkinje) {
 
     struct ensight_grid *ensight_grid;
 
     if(grid->num_active_cells > 0 && grid->purkinje) {
 
-        uint32_t num_active_cells = grid->num_active_cells;
-        uint32_t number_of_purkinje_cells = grid->purkinje->num_active_purkinje_cells;
-        ensight_grid = new_ensight_grid(2);
+        if (save_purkinje) {
+            uint32_t num_active_cells = grid->num_active_cells;
+            uint32_t number_of_purkinje_cells = grid->purkinje->num_active_purkinje_cells;
+            ensight_grid = new_ensight_grid(2);
 
-        arrsetcap(ensight_grid->parts[0].points, num_active_cells * 8);
-        arrsetcap(ensight_grid->parts[0].cell_visibility, num_active_cells);
-        arrsetcap(ensight_grid->parts[0].cells,  num_active_cells);
+            arrsetcap(ensight_grid->parts[0].points, num_active_cells * 8);
+            arrsetcap(ensight_grid->parts[0].cell_visibility, num_active_cells);
+            arrsetcap(ensight_grid->parts[0].cells,  num_active_cells);
 
-        arrsetcap(ensight_grid->parts[1].points, number_of_purkinje_cells * 2);
-        arrsetcap(ensight_grid->parts[1].cells,  number_of_purkinje_cells);
+            arrsetcap(ensight_grid->parts[1].points, number_of_purkinje_cells * 2);
+            arrsetcap(ensight_grid->parts[1].cells,  number_of_purkinje_cells);
 
-        ensight_grid->max_v = FLT_MIN;
-        ensight_grid->min_v = FLT_MAX;
+            ensight_grid->max_v = FLT_MIN;
+            ensight_grid->min_v = FLT_MAX;
+        }
+        else {
+            ensight_grid = new_ensight_grid(1);
 
+            uint32_t num_active_cells = grid->num_active_cells;
+            arrsetcap(ensight_grid->parts[0].points, num_active_cells * 8);
+            arrsetcap(ensight_grid->parts[0].cell_visibility, num_active_cells);
+            arrsetcap(ensight_grid->parts[0].cells,  num_active_cells);
+        }
     } else {
         ensight_grid = new_ensight_grid(1);
 
@@ -399,7 +408,6 @@ struct ensight_grid * new_ensight_grid_from_alg_grid(struct grid *grid, bool cli
             arrsetcap(ensight_grid->parts[0].points, number_of_purkinje_cells * 2);
             arrsetcap(ensight_grid->parts[0].cells,  number_of_purkinje_cells);
         }
-
     }
 
     ensight_grid->max_v = FLT_MIN;
@@ -511,7 +519,7 @@ struct ensight_grid * new_ensight_grid_from_alg_grid(struct grid *grid, bool cli
         hash = NULL;
     }
 
-    if(grid->purkinje) {
+    if(grid->purkinje && save_purkinje) {
 
         int part_n = 0;
 

src/ensight_utils/ensight_grid.h

@@ -33,12 +33,12 @@ struct ensight_grid * new_ensight_grid(uint32_t num_parts);
 struct ensight_grid * new_ensight_grid_from_alg_grid(struct grid *grid, bool clip_with_plain,
                                                                      float *plain_coordinates, bool clip_with_bounds,
                                                                      float *bounds, bool read_fibers_f,
-                                                                     bool save_fibrotic);
+                                                                     bool save_fibrotic, bool save_purkinje);
 
 void save_ensight_grid_as_ensight6_geometry(struct ensight_grid *grid, char *filename, bool binary);
 
 void free_ensight_grid(struct ensight_grid *ensight_grid);
 void save_case_file(char *filename, uint64_t num_files, real_cpu dt, int print_rate, int num_state_var);
-void save_en6_result_file(char *filename, struct grid *the_grid, bool binary);
+void save_en6_result_file(char *filename, struct grid *the_grid, bool binary, bool save_purkinje);
 void save_en6_result_file_state_vars(char *filename, real *sv_cpu, size_t num_cells, size_t num_sv_entries, int sv_entry, bool binary, bool gpu);
 #endif //MONOALG3D_ENSIGHT_GRID_H