Add missing files

src/matrix_assembly_library/matrix_assembly.c

@@ -1220,3 +1220,79 @@ 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);
+}