From 3abede798ed8ffed5762f2177d2d4511c1363780 Mon Sep 17 00:00:00 2001 From: bergolho Date: Tue, 5 Nov 2024 10:58:16 +0000 Subject: [PATCH] Add 'custom_functions' files --- src/domains_library/custom_functions.c | 136 ++++++ src/domains_library/custom_mesh_info_data.h | 9 + src/extra_data_library/custom_functions.c | 40 ++ .../custom_functions.c | 433 ++++++++++++++++++ 4 files changed, 618 insertions(+) create mode 100644 src/domains_library/custom_functions.c create mode 100644 src/domains_library/custom_mesh_info_data.h create mode 100644 src/extra_data_library/custom_functions.c create mode 100644 src/matrix_assembly_library/custom_functions.c diff --git a/src/domains_library/custom_functions.c b/src/domains_library/custom_functions.c new file mode 100644 index 00000000..2b35a6b3 --- /dev/null +++ b/src/domains_library/custom_functions.c @@ -0,0 +1,136 @@ +// +// Created by bergolho and jennyhelyanwe on 03/06/24. +// + +#include "domain_helpers.h" + +#include "../3dparty/sds/sds.h" +#include "../3dparty/stb_ds.h" +#include "../config/domain_config.h" +#include "../config_helpers/config_helpers.h" +#include "../libraries_common/common_data_structures.h" +#include "../logger/logger.h" +#include "../utils/stop_watch.h" +#include "../utils/utils.h" +#include +#include + +#include "mesh_info_data.h" +#include "custom_mesh_info_data.h" + +#include +#include + +SET_CUSTOM_DATA_FOR_MESH(set_custom_data_for_dti_mesh_with_dense_and_sparse_regions_twave) { + + INITIALIZE_DTI_MESH_INFO(cell); + + // Ruben mesh labels (transmuralityLabels): 1=ENDO, 2=MID, 3=EPI + // Dense/Sparse regions (DenseSparseRegion): 1=NORMAL, 2=DENSE, 3=SPARSE + // Cellular model labels: 0=ENDO, 1=MID, 2=EPI, 3=DENSE, 4=SPARSE + DTI_MESH_TRANSMURALITY(cell) = custom_data[0]; + DTI_MESH_TRANSMURALITY_LABELS(cell) = custom_data[1]-1; + DTI_MESH_BASE_APEX_HETEROGENEITY(cell) = custom_data[2]; + DTI_MESH_FAST_ENDO(cell) = (int)custom_data[3]; + + cell->sigma.fibers.f[0] = custom_data[5]; + cell->sigma.fibers.f[1] = custom_data[6]; + cell->sigma.fibers.f[2] = custom_data[7]; + + cell->sigma.fibers.s[0] = custom_data[8]; + cell->sigma.fibers.s[1] = custom_data[9]; + cell->sigma.fibers.s[2] = custom_data[10]; + + cell->sigma.fibers.n[0] = custom_data[11]; + cell->sigma.fibers.n[1] = custom_data[12]; + cell->sigma.fibers.n[2] = custom_data[13]; + + DTI_MESH_SCALE_FACTOR_IKS(cell) = custom_data[14]; +} + +// Domain function to load the ALG file +// Works with all DTI meshes (T-wave personalisation) - Fiber orientation, transmurality, dense/sparse regions +// and scale factor for IKs (resolution = 500um) +SET_SPATIAL_DOMAIN(initialize_grid_with_dti_mesh_with_dense_and_sparse_regions_twave) { + + char *mesh_file = NULL; + GET_PARAMETER_STRING_VALUE_OR_REPORT_ERROR(mesh_file, config, "mesh_file"); + + real_cpu start_h = 500.0; + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real_cpu, start_h, config, "original_discretization"); + + real_cpu desired_h = 500.0; + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real_cpu, desired_h, config, "desired_discretization"); + + assert(the_grid); + + log_info("DTI Mesh with discretization: %lf\n", desired_h); + + uint32_t num_volumes = 0; + GET_PARAMETER_NUMERIC_VALUE_OR_REPORT_ERROR(uint32_t, num_volumes, config, "num_volumes"); + + // ALG mesh column config: + // center_x, center_y, center_z, dx, dy, dz, transmurality, transmurality_label, base_to_apex, dense_sparse_region_tag, fibrotic, fx, fy, fz, sx, sy, sz, nx, ny, nz, sf_IKs + uint32_t num_extra_fields = 15; + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(uint32_t, num_extra_fields, config, "num_extra_fields"); + + int num_loaded = set_custom_mesh_from_file(the_grid, mesh_file, num_volumes, start_h, num_extra_fields, set_custom_data_for_dti_mesh_with_dense_and_sparse_regions_twave); + + log_info("Read %d volumes from file (expected %d): %s\n", num_loaded, num_volumes, mesh_file); + + int num_refs_remaining = (int)(start_h / desired_h) - 1; + + if(num_refs_remaining > 0) { + refine_grid(the_grid, num_refs_remaining); + } + + the_grid->start_discretization = SAME_POINT3D(desired_h); + + free(mesh_file); + + return num_loaded; +} + +SET_SPATIAL_DOMAIN(initialize_grid_with_dti_mesh_with_dense_and_sparse_regions_twave_mesh_number) { + + char *mesh_directory = NULL; + GET_PARAMETER_STRING_VALUE_OR_REPORT_ERROR(mesh_directory, config, "mesh_directory"); + int mesh_number = 0; + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(int, mesh_number, config, "mesh_number"); + + real_cpu start_h = 500.0; + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real_cpu, start_h, config, "original_discretization"); + + real_cpu desired_h = 500.0; + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real_cpu, desired_h, config, "desired_discretization"); + + assert(the_grid); + + log_info("DTI Mesh with discretization: %lf\n", desired_h); + + char mesh_file[500]; + sprintf(mesh_file,"%s/%d.alg", mesh_directory, mesh_number); + uint32_t num_volumes = 0; + GET_PARAMETER_NUMERIC_VALUE_OR_REPORT_ERROR(uint32_t, num_volumes, config, "num_volumes"); + + // ALG mesh column config: + // center_x, center_y, center_z, dx, dy, dz, transmurality, transmurality_label, base_to_apex, dense_sparse_region_tag, fibrotic, fx, fy, fz, sx, sy, sz, nx, ny, nz, sf_IKs + uint32_t num_extra_fields = 15; + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(uint32_t, num_extra_fields, config, "num_extra_fields"); + + int num_loaded = set_custom_mesh_from_file(the_grid, mesh_file, num_volumes, start_h, num_extra_fields, set_custom_data_for_dti_mesh_with_dense_and_sparse_regions_twave); + + log_info("Read %d volumes from file (expected %d): %s\n", num_loaded, num_volumes, mesh_file); + + int num_refs_remaining = (int)(start_h / desired_h) - 1; + + if(num_refs_remaining > 0) { + refine_grid(the_grid, num_refs_remaining); + } + + the_grid->start_discretization = SAME_POINT3D(desired_h); + + free(mesh_file); + + return num_loaded; +} diff --git a/src/domains_library/custom_mesh_info_data.h b/src/domains_library/custom_mesh_info_data.h new file mode 100644 index 00000000..7b18f1b8 --- /dev/null +++ b/src/domains_library/custom_mesh_info_data.h @@ -0,0 +1,9 @@ +#ifndef __CUSTOM_MESH_INFO_DATA_H +#define __CUSTOM_MESH_INFO_DATA_H + +#include +#include "../common_types/common_types.h" + + + +#endif /* __CUSTOM_MESH_INFO_DATA_H */ diff --git a/src/extra_data_library/custom_functions.c b/src/extra_data_library/custom_functions.c new file mode 100644 index 00000000..9015cb6e --- /dev/null +++ b/src/extra_data_library/custom_functions.c @@ -0,0 +1,40 @@ +// +// Created by bergolho and jennyhelyanwe on 03/06/24. +// + +#include + +#include "../config/extra_data_config.h" +#include "../config_helpers/config_helpers.h" +#include "../libraries_common/common_data_structures.h" +#include "../domains_library/mesh_info_data.h" +#include "../domains_library/custom_mesh_info_data.h" +#include "../logger/logger.h" +#include "helper_functions.h" + +// Extra data function to configure the transmurality and IKs scaling factor +// Works with all DTI meshes - Extra data function to load the ToRORd cellular model with modifications +// Function used for the T-wave personalisation paper. +SET_EXTRA_DATA(set_extra_data_for_dti_mesh_twave_with_torord_gksgkrtjca) { + + uint32_t num_active_cells = the_grid->num_active_cells; + struct cell_node ** ac = the_grid->active_cells; + + struct extra_data_for_torord_gksgkrtjca_twave *extra_data = NULL; + extra_data = set_common_torord_gksgkrtjca_twave_data(config, num_active_cells); // ToRORd used for T-wave personalisation + + // Ruben mesh labels (transmuralityLabels): 1=ENDO, 2=MID, 3=EPI + // Dense/Sparse regions (DenseSparseRegion): 1=NORMAL, 2=DENSE, 3=SPARSE + // Cellular model labels: 0=ENDO, 1=MID, 2=EPI, 3=DENSE, 4=SPARSE + OMP(parallel for) + for (int i = 0; i < num_active_cells; i++) { + real tag = DTI_MESH_TRANSMURALITY_LABELS(ac[i]); + real sf_IKs = DTI_MESH_SCALE_FACTOR_IKS(ac[i]); + extra_data->transmurality[i] = tag; + extra_data->sf_IKs[i] = sf_IKs; + } + + SET_EXTRA_DATA_SIZE(sizeof(struct extra_data_for_torord_gksgkrtjca_twave)); + + return (void*)extra_data; +} diff --git a/src/matrix_assembly_library/custom_functions.c b/src/matrix_assembly_library/custom_functions.c new file mode 100644 index 00000000..29f5538d --- /dev/null +++ b/src/matrix_assembly_library/custom_functions.c @@ -0,0 +1,433 @@ +// +// Created by bergolho and jennyhelyanwe on 03/06/24. +// + +#include +#include +#include +#include +#include + +#include "../3dparty/sds/sds.h" +#include "../3dparty/stb_ds.h" +#include "../alg/grid/grid.h" +#include "../config/assembly_matrix_config.h" +#include "../libraries_common/common_data_structures.h" +#include "../utils/file_utils.h" +#include "../utils/utils.h" +#include "../domains_library/custom_mesh_info_data.h" +#include "purkinje_coupling_matrix_assembly.c" + +#ifdef COMPILE_CUDA +#include "../gpu_utils/gpu_utils.h" +#endif + +// Matrix assembly function to configure the conductivities for the NORMAL, DENSE and SPARSE regions. +// Works with all DTI meshes +// Function used for the T-wave personalisation paper. +ASSEMBLY_MATRIX(anisotropic_sigma_assembly_matrix_with_dense_and_sparse_region_twave) { + + if(the_grid->adaptive) { + log_error_and_exit("anisotropic_sigma_assembly_matrix_with_dense_and_sparse_region_twave function does not support mesh adaptivity yet!. Aborting!\n"); + } + + 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); + + // D tensor // + // | sx sxy sxz | + // | sxy sy syz | + // | sxz syz sz | + real_cpu D[3][3]; + real_cpu D_dense[3][3]; + real_cpu D_sparse[3][3]; + int i; + + real_cpu sigma_l = 0.0; + real_cpu sigma_t = 0.0; + real_cpu sigma_n = 0.0; + real_cpu sigma_l_dense = 0.0; + real_cpu sigma_t_dense = 0.0; + real_cpu sigma_n_dense = 0.0; + real_cpu sigma_l_sparse = 0.0; + real_cpu sigma_t_sparse = 0.0; + real_cpu sigma_n_sparse = 0.0; + + char *fiber_file = NULL; + GET_PARAMETER_STRING_VALUE_OR_USE_DEFAULT(fiber_file, config, "fibers_file"); + + bool fibers_in_mesh = false; + GET_PARAMETER_BOOLEAN_VALUE_OR_USE_DEFAULT(fibers_in_mesh, config, "fibers_in_mesh"); + + struct fiber_coords *fibers = NULL; + + // Normal myocardium conductivities + 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"); + + // Dense region myocardium conductivities + GET_PARAMETER_NUMERIC_VALUE_OR_REPORT_ERROR(real_cpu, sigma_l_dense, config, "sigma_l_dense"); + GET_PARAMETER_NUMERIC_VALUE_OR_REPORT_ERROR(real_cpu, sigma_t_dense, config, "sigma_t_dense"); + GET_PARAMETER_NUMERIC_VALUE_OR_REPORT_ERROR(real_cpu, sigma_n_dense, config, "sigma_n_dense"); + + // Sparse myocardium conductivities + GET_PARAMETER_NUMERIC_VALUE_OR_REPORT_ERROR(real_cpu, sigma_l_sparse, config, "sigma_l_sparse"); + GET_PARAMETER_NUMERIC_VALUE_OR_REPORT_ERROR(real_cpu, sigma_t_sparse, config, "sigma_t_sparse"); + GET_PARAMETER_NUMERIC_VALUE_OR_REPORT_ERROR(real_cpu, sigma_n_sparse, config, "sigma_n_sparse"); + + real_cpu *f = NULL; + real_cpu *s = NULL; + real_cpu *n = NULL; + + if(fiber_file) { + log_info("Loading mesh fibers\n"); + fibers = read_fibers(fiber_file, false); + } + else if(!fibers_in_mesh) { + GET_PARAMETER_VECTOR_VALUE_OR_USE_DEFAULT(f, config, "f", 3); + GET_PARAMETER_VECTOR_VALUE_OR_USE_DEFAULT(s, config, "s", 3); + GET_PARAMETER_VECTOR_VALUE_OR_USE_DEFAULT(n, config, "n", 3); + + if(!f) { + f = malloc(sizeof(real_cpu)*3); + f[0] = 1.0; + f[1] = 0.0; + f[2] = 0.0; + } + + if(!s) { + s = malloc(sizeof(real_cpu)*3); + s[0] = 0.0; + s[1] = 1.0; + s[2] = 0.0; + } + + if(!n) { + n = malloc(sizeof(real_cpu)*3); + n[0] = 0.0; + n[1] = 0.0; + n[2] = 1.0; + } + } + + OMP(parallel for private(D)) + for(i = 0; i < num_active_cells; i++) { + + if(fibers) { + 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, sigma_t); + calc_tensor2(D_dense, fibers[fiber_index].f, sigma_l_dense, sigma_t_dense); + calc_tensor2(D_sparse, fibers[fiber_index].f, sigma_l_sparse, sigma_t_sparse); + } + else { + calc_tensor(D, fibers[fiber_index].f, fibers[fiber_index].s, fibers[fiber_index].n, sigma_l, sigma_t, sigma_n); + calc_tensor(D_dense, fibers[fiber_index].f, fibers[fiber_index].s, fibers[fiber_index].n, sigma_l_dense, sigma_t_dense, sigma_n_dense); + calc_tensor(D_sparse, fibers[fiber_index].f, fibers[fiber_index].s, fibers[fiber_index].n, sigma_l_sparse, sigma_t_sparse, sigma_n_sparse); + } + ac[i]->sigma.fibers = fibers[fiber_index]; + } + else if(fibers_in_mesh) { + if(sigma_t == sigma_n) { + calc_tensor2(D, ac[i]->sigma.fibers.f, sigma_l, sigma_t); + calc_tensor2(D_dense, ac[i]->sigma.fibers.f, sigma_l_dense, sigma_t_dense); + calc_tensor2(D_sparse, ac[i]->sigma.fibers.f, sigma_l_sparse, sigma_t_sparse); + } + else { + calc_tensor(D, ac[i]->sigma.fibers.f, ac[i]->sigma.fibers.s, ac[i]->sigma.fibers.n, sigma_l, sigma_t, sigma_n); + calc_tensor(D_dense, ac[i]->sigma.fibers.f, ac[i]->sigma.fibers.s, ac[i]->sigma.fibers.n, sigma_l_dense, sigma_t_dense, sigma_n_dense); + calc_tensor(D_sparse, ac[i]->sigma.fibers.f, ac[i]->sigma.fibers.s, ac[i]->sigma.fibers.n, sigma_l_sparse, sigma_t_sparse, sigma_n_sparse); + } + + } + else { + if(sigma_t == sigma_n) { + calc_tensor2(D, f, sigma_l, sigma_t); + calc_tensor2(D_dense, f, sigma_l_dense, sigma_t_dense); + calc_tensor2(D_sparse, f, sigma_l_sparse, sigma_t_sparse); + } + else { + calc_tensor(D, f, s, n, sigma_l, sigma_t, sigma_n); + calc_tensor(D_dense, f, s, n, sigma_l_dense, sigma_t_dense, sigma_n_dense); + calc_tensor(D_sparse, f, s, n, sigma_l_sparse, sigma_t_sparse, sigma_n_sparse); + } + } + + // Check if the cell is tagged with an endocardium region tag + // 2 = DENSE region || 3 = SPARSE region + int tag = DTI_MESH_FAST_ENDO(ac[i]); + + // Cell is inside the DENSE endocardium region + if (tag == 2) { + ac[i]->sigma.x = D_dense[0][0]; + ac[i]->sigma.y = D_dense[1][1]; + ac[i]->sigma.z = D_dense[2][2]; + + ac[i]->sigma.xy = D_dense[0][1]; + ac[i]->sigma.xz = D_dense[0][2]; + ac[i]->sigma.yz = D_dense[1][2]; + } + // Cell is inside the SPARSE endocardium region + else if (tag == 3) { + ac[i]->sigma.x = D_sparse[0][0]; + ac[i]->sigma.y = D_sparse[1][1]; + ac[i]->sigma.z = D_sparse[2][2]; + + ac[i]->sigma.xy = D_sparse[0][1]; + ac[i]->sigma.xz = D_sparse[0][2]; + ac[i]->sigma.yz = D_sparse[1][2]; + } + // Otherwise, NORMAL type of cell + else { + 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); +} + +// Matrix assembly function to configure the conductivities for the NORMAL, DENSE and SPARSE regions. +// + Works with Purkinje coupled simulations. +// Works with all DTI meshes +// Function used for the T-wave personalisation paper. +ASSEMBLY_MATRIX(purkinje_coupling_with_anisotropic_sigma_assembly_matrix_with_dense_and_sparse_region_twave) { + + if(the_grid->adaptive) { + log_error_and_exit("purkinje_coupling_with_anisotropic_sigma_assembly_matrix_with_dense_and_sparse_region_twave function does not support mesh adaptivity yet!. Aborting!\n"); + } + +// [MYOCARDIUM SETUP] + 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); + + // D tensor // + // | sx sxy sxz | + // | sxy sy syz | + // | sxz syz sz | + real_cpu D[3][3]; + real_cpu D_dense[3][3]; + real_cpu D_sparse[3][3]; + int i; + + real_cpu sigma_l = 0.0; + real_cpu sigma_t = 0.0; + real_cpu sigma_n = 0.0; + real_cpu sigma_l_dense = 0.0; + real_cpu sigma_t_dense = 0.0; + real_cpu sigma_n_dense = 0.0; + real_cpu sigma_l_sparse = 0.0; + real_cpu sigma_t_sparse = 0.0; + real_cpu sigma_n_sparse = 0.0; + real_cpu sigma_purkinje = 0.0; + + char *fiber_file = NULL; + GET_PARAMETER_STRING_VALUE_OR_USE_DEFAULT(fiber_file, config, "fibers_file"); + + bool fibers_in_mesh = false; + GET_PARAMETER_BOOLEAN_VALUE_OR_USE_DEFAULT(fibers_in_mesh, config, "fibers_in_mesh"); + + struct fiber_coords *fibers = NULL; + + // Normal myocardium conductivities + 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"); + + // Dense region myocardium conductivities + GET_PARAMETER_NUMERIC_VALUE_OR_REPORT_ERROR(real_cpu, sigma_l_dense, config, "sigma_l_dense"); + GET_PARAMETER_NUMERIC_VALUE_OR_REPORT_ERROR(real_cpu, sigma_t_dense, config, "sigma_t_dense"); + GET_PARAMETER_NUMERIC_VALUE_OR_REPORT_ERROR(real_cpu, sigma_n_dense, config, "sigma_n_dense"); + + // Sparse myocardium conductivities + GET_PARAMETER_NUMERIC_VALUE_OR_REPORT_ERROR(real_cpu, sigma_l_sparse, config, "sigma_l_sparse"); + GET_PARAMETER_NUMERIC_VALUE_OR_REPORT_ERROR(real_cpu, sigma_t_sparse, config, "sigma_t_sparse"); + GET_PARAMETER_NUMERIC_VALUE_OR_REPORT_ERROR(real_cpu, sigma_n_sparse, config, "sigma_n_sparse"); + + // Purkinje conductivity + GET_PARAMETER_NUMERIC_VALUE_OR_USE_DEFAULT(real, sigma_purkinje, config, "sigma_purkinje"); + + real_cpu *f = NULL; + real_cpu *s = NULL; + real_cpu *n = NULL; + + if(fiber_file) { + log_info("Loading mesh fibers\n"); + fibers = read_fibers(fiber_file, false); + } + else if(!fibers_in_mesh) { + GET_PARAMETER_VECTOR_VALUE_OR_USE_DEFAULT(f, config, "f", 3); + GET_PARAMETER_VECTOR_VALUE_OR_USE_DEFAULT(s, config, "s", 3); + GET_PARAMETER_VECTOR_VALUE_OR_USE_DEFAULT(n, config, "n", 3); + + if(!f) { + f = malloc(sizeof(real_cpu)*3); + f[0] = 1.0; + f[1] = 0.0; + f[2] = 0.0; + } + + if(!s) { + s = malloc(sizeof(real_cpu)*3); + s[0] = 0.0; + s[1] = 1.0; + s[2] = 0.0; + } + + if(!n) { + n = malloc(sizeof(real_cpu)*3); + n[0] = 0.0; + n[1] = 0.0; + n[2] = 1.0; + } + } + + OMP(parallel for private(D)) + for(i = 0; i < num_active_cells; i++) { + + if(fibers) { + 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, sigma_t); + calc_tensor2(D_dense, fibers[fiber_index].f, sigma_l_dense, sigma_t_dense); + calc_tensor2(D_sparse, fibers[fiber_index].f, sigma_l_sparse, sigma_t_sparse); + } + else { + calc_tensor(D, fibers[fiber_index].f, fibers[fiber_index].s, fibers[fiber_index].n, sigma_l, sigma_t, sigma_n); + calc_tensor(D_dense, fibers[fiber_index].f, fibers[fiber_index].s, fibers[fiber_index].n, sigma_l_dense, sigma_t_dense, sigma_n_dense); + calc_tensor(D_sparse, fibers[fiber_index].f, fibers[fiber_index].s, fibers[fiber_index].n, sigma_l_sparse, sigma_t_sparse, sigma_n_sparse); + } + ac[i]->sigma.fibers = fibers[fiber_index]; + } + else if(fibers_in_mesh) { + if(sigma_t == sigma_n) { + calc_tensor2(D, ac[i]->sigma.fibers.f, sigma_l, sigma_t); + calc_tensor2(D_dense, ac[i]->sigma.fibers.f, sigma_l_dense, sigma_t_dense); + calc_tensor2(D_sparse, ac[i]->sigma.fibers.f, sigma_l_sparse, sigma_t_sparse); + } + else { + calc_tensor(D, ac[i]->sigma.fibers.f, ac[i]->sigma.fibers.s, ac[i]->sigma.fibers.n, sigma_l, sigma_t, sigma_n); + calc_tensor(D_dense, ac[i]->sigma.fibers.f, ac[i]->sigma.fibers.s, ac[i]->sigma.fibers.n, sigma_l_dense, sigma_t_dense, sigma_n_dense); + calc_tensor(D_sparse, ac[i]->sigma.fibers.f, ac[i]->sigma.fibers.s, ac[i]->sigma.fibers.n, sigma_l_sparse, sigma_t_sparse, sigma_n_sparse); + } + + } + else { + if(sigma_t == sigma_n) { + calc_tensor2(D, f, sigma_l, sigma_t); + calc_tensor2(D_dense, f, sigma_l_dense, sigma_t_dense); + calc_tensor2(D_sparse, f, sigma_l_sparse, sigma_t_sparse); + } + else { + calc_tensor(D, f, s, n, sigma_l, sigma_t, sigma_n); + calc_tensor(D_dense, f, s, n, sigma_l_dense, sigma_t_dense, sigma_n_dense); + calc_tensor(D_sparse, f, s, n, sigma_l_sparse, sigma_t_sparse, sigma_n_sparse); + } + } + + // Check if the cell is tagged with an endocardium region tag + // 2 = DENSE region || 3 = SPARSE region + int tag = DTI_MESH_FAST_ENDO(ac[i]); + + // Cell is inside the DENSE endocardium region + if (tag == 2) { + ac[i]->sigma.x = D_dense[0][0]; + ac[i]->sigma.y = D_dense[1][1]; + ac[i]->sigma.z = D_dense[2][2]; + + ac[i]->sigma.xy = D_dense[0][1]; + ac[i]->sigma.xz = D_dense[0][2]; + ac[i]->sigma.yz = D_dense[1][2]; + } + // Cell is inside the SPARSE endocardium region + else if (tag == 3) { + ac[i]->sigma.x = D_sparse[0][0]; + ac[i]->sigma.y = D_sparse[1][1]; + ac[i]->sigma.z = D_sparse[2][2]; + + ac[i]->sigma.xy = D_sparse[0][1]; + ac[i]->sigma.xz = D_sparse[0][2]; + ac[i]->sigma.yz = D_sparse[1][2]; + } + // Otherwise, NORMAL type of cell + else { + 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); + +// [PURKINJE SETUP] + static bool sigma_purkinje_initialized = false; + + uint32_t num_purkinje_active_cells = the_grid->purkinje->num_active_purkinje_cells; + struct cell_node **ac_purkinje = the_grid->purkinje->purkinje_cells; + struct node *pk_node = the_grid->purkinje->network->list_nodes; + bool has_point_data = the_grid->purkinje->network->has_point_data; + + initialize_diagonal_elements_purkinje(the_solver, the_grid); + + if(!sigma_purkinje_initialized) { + // Check if the Purkinje network file has the POINT_DATA section + if (has_point_data) { + struct node *tmp = the_grid->purkinje->network->list_nodes; + uint32_t i = 0; + while (tmp != NULL) + { + // Copy the prescribed conductivity from the Purkinje network file into the ALG Purkinje cell structure + ac_purkinje[i]->sigma.x = tmp->sigma; + + tmp = tmp->next; i++; + } + } + // Otherwise, initilize the conductivity of all cells homogenously with the value from the configuration file + else { + OMP(parallel for) + for (uint32_t i = 0; i < num_active_cells; i++) { + ac[i]->sigma.x = sigma_purkinje; + } + } + sigma_purkinje_initialized = true; + } + fill_discretization_matrix_elements_purkinje(has_point_data,sigma_purkinje,ac_purkinje,num_purkinje_active_cells,pk_node); +}