adding fiber generation codes

utilities/fiber_generation/README.md

@@ -0,0 +1,35 @@
+# SV-fibergen
+Python + SVmultiphysics codes for fiber generation. Two methods are implemented:
+* Bayer et al. (2012). [link](https://doi.org/10.1007/s10439-012-0593-5)
+* Doste et al. (2018). [link](https://doi.org/10.1002/cnm.3185)
+
+### Examples
+The `main_bayer.py` and `main_doste.py` are scripts to run both methods in the geometry described in the `example/truncated` and `example/ot` folders respectively.
+
+<img src="example/truncated/bayer.png" alt="Results for truncated BiV (Bayer)" width="640" />
+<img src="example/ot/doste.png" alt="Results for BiV w/ outflow tracts (Doste)" width="640" />
+
+Note that the Doste methods needs a geometry with outflow tracts to be run (each valve needs to be defined as a separated surface). Bayer can be run in any biventricular geometry.
+
+
+### Updates to the old code
+* All operations are vectorized now.
+* The SVmultiphysics solver now solves a Laplace equation.
+* In Bayer: For the bislerp interpolation, instead of using the correction described in Bayer et al. (that returns a discontinuity), the basis are flipped to maintain a coherent fiber direction (see function `generate_fibers_BiV_Bayer_cells` in `FibGen.py`).
+* In Bayer: The beta angles were not being included correctly. The second rotation was being applied respect the first vector (circumferential) when it should be respect to the second vector (longitudinal) (see function `generate_fibers_BiV_Bayer_cells` in `FibGen.py`).
+
+### Notes on SVmultiphysics solver
+
+To solve a Laplace equation directly from the transient HEAT solver in SVmultiphysics, in `<GeneralSimulationParameters>` we need to set,
+```
+<Number_of_time_steps> 1 </Number_of_time_steps>
+<Time_step_size> 1 </Time_step_size>
+<Spectral_radius_of_infinite_time_step> 0. </Spectral_radius_of_infinite_time_step>
+```
+and in `<Add_equation type="heatS" >`,
+```
+<Conductivity> 1.0 </Conductivity>
+<Source_term> 0.0 </Source_term>
+<Density> 0.0 </Density>
+```
+This will allow us to solve the Laplace equation directly in 1 timestep and 1 iteration.

utilities/fiber_generation/main_bayer.py

@@ -0,0 +1,64 @@
+#!/usr/bin/env python
+# -*-coding:utf-8 -*-
+'''
+Created on 2025/11/21 20:38:14
+
+@author: Javiera Jilberto Vallejos 
+'''
+
+import os
+import src.FibGen as fg
+from time import time
+
+###########################################################
+############  USER INPUTS  ################################
+###########################################################
+
+run_flag = True
+svfsi_exec = "svmultiphysics "
+
+mesh_path = "example/truncated/VOLUME.vtu"
+surfaces_dir = f"example/truncated/mesh-surfaces"
+outdir = "example/truncated/output_b"
+
+surface_names = {'epi': 'EPI.vtp',
+                 'epi_apex': 'EPI_APEX.vtp',    # New surface
+                 'base': 'BASE.vtp',
+                 'endo_lv': 'LV.vtp',
+                 'endo_rv': 'RV.vtp'}
+
+# Parameters for the Bayer et al. method https://doi.org/10.1007/s10439-012-0593-5. 
+params = {
+    "ALFA_END": 60.0,
+    "ALFA_EPI": -60.0,
+    "BETA_END": 20.0,
+    "BETA_EPI": -20.0,
+}
+
+
+###########################################################
+############  FIBER GENERATION  ###########################
+###########################################################
+
+# Make sure the paths are full paths
+mesh_path = os.path.abspath(mesh_path)
+surfaces_dir = os.path.abspath(surfaces_dir)
+outdir = os.path.abspath(outdir)
+
+start = time()
+fg.generate_epi_apex(mesh_path, surfaces_dir, surface_names)
+
+# Run the Laplace solver
+if run_flag:
+    template_file = os.path.join(os.path.dirname(os.path.abspath(__file__)), "src", "templates", "solver_bayer.xml")
+    laplace_results_file = fg.runLaplaceSolver(mesh_path, surfaces_dir, mesh_path, svfsi_exec, template_file, outdir, surface_names)
+laplace_results_file = outdir + '/result_001.vtu'
+
+# Generate the fiber directions
+result_mesh = fg.generate_fibers_BiV_Bayer_cells(outdir, laplace_results_file, params, return_angles=True, return_intermediate=True)
+
+print(f"generate fibers (Bayer method) elapsed time: {time() - start:.3f} s")
+
+# Optional, save the result mesh with intermediate field and angles for checking
+result_mesh_path = os.path.join(outdir, "results_bayer.vtu")
+result_mesh.save(result_mesh_path)

utilities/fiber_generation/main_doste.py

@@ -0,0 +1,85 @@
+#!/usr/bin/env python
+# -*-coding:utf-8 -*-
+'''
+Created on 2025/11/21 20:38:14
+
+@author: Javiera Jilberto Vallejos 
+'''
+
+import os
+import src.FibGen as fg
+from time import time
+
+###########################################################
+############  USER INPUTS  ################################
+###########################################################
+
+run_flag = True
+method = 'doste'
+svfsi_exec = "svmultiphysics "
+
+mesh_path = "example/ot/mesh-complete.mesh.vtu"
+surfaces_dir = f"example/ot/mesh-surfaces"
+outdir = "example/ot/output_d"
+
+surface_names = {'epi': 'epi.vtp',
+                 'epi_apex': 'epi_apex.vtp',    # New surface
+                 'av': 'av.vtp',
+                 'mv': 'mv.vtp',
+                 'tv': 'tv.vtp',
+                 'pv': 'pv.vtp',
+                 'base': 'top.vtp',             # This is all the valves together, it is used to find the apex.
+                 'endo_lv': 'endo_lv.vtp',
+                 'endo_rv': 'endo_rv.vtp'}
+
+# Parameters from the Doste paper https://doi.org/10.1002/cnm.3185
+params = {
+    # A = alpha angle
+    'AENDORV' : 90,
+    'AEPIRV' : -25,
+    'AENDOLV' : 60,
+    'AEPILV' : -60,
+
+    'AOTENDOLV' : 90, 
+    'AOTENDORV' : 90,
+    'AOTEPILV' : 0,
+    'AOTEPIRV' : 0,
+
+    # B = beta angle (this have an opposite sign to the Doste paper, 
+    # but it's because the longitudinal direction is opposite)
+    'BENDORV' : 20,
+    'BEPIRV' : -20,
+    'BENDOLV' : 20,
+    'BEPILV' : -20,
+}
+
+
+###########################################################
+############  FIBER GENERATION  ###########################
+###########################################################
+
+# Make sure the paths are full paths
+mesh_path = os.path.abspath(mesh_path)
+surfaces_dir = os.path.abspath(surfaces_dir)
+outdir = os.path.abspath(outdir)
+
+# Generate the apex surface
+start = time()
+
+start = time()
+fg.generate_epi_apex(mesh_path, surfaces_dir, surface_names)
+
+# Run the Laplace solver
+if run_flag:
+    template_file = os.path.join(os.path.dirname(os.path.abspath(__file__)), "src", "templates", "solver_doste.xml")
+    laplace_results_file = fg.runLaplaceSolver(mesh_path, surfaces_dir, mesh_path, svfsi_exec, template_file, outdir, surface_names)
+laplace_results_file = outdir + '/result_001.vtu'
+
+# Generate the fiber directions
+result_mesh = fg.generate_fibers_BiV_Doste_cells(outdir, laplace_results_file, params, return_angles=True, return_intermediate=False)
+
+print(f"generate fibers (Doste method) elapsed time: {time() - start:.3f} s")
+
+# Optional, save the result mesh with intermediate field and angles for checking
+result_mesh_path = os.path.join(outdir, "results_doste.vtu")
+result_mesh.save(result_mesh_path)