Add VTU to XDMF mesh conversion using meshio

environment.yml

@@ -12,6 +12,7 @@ dependencies:
   - pip
   - python >=3.8
   - git
+  - meshio
   - pip:
       - git+https://github.com/KVSlab/VaMPy.git
       - git+https://github.com/KVSlab/turtleFSI.git

src/fsipy/automatedPreprocessing/automated_preprocessing.py

@@ -14,13 +14,13 @@
 from vampy.automatedPreprocessing import ToolRepairSTL
 from vampy.automatedPreprocessing.preprocessing_common import read_polydata, get_centers_for_meshing, \
     dist_sphere_diam, dist_sphere_curvature, dist_sphere_constant, get_regions_to_refine, add_flow_extension, \
-    write_mesh, mesh_alternative, find_boundaries, compute_flow_rate, setup_model_network, \
+    mesh_alternative, find_boundaries, compute_flow_rate, setup_model_network, \
     radiusArrayName, scale_surface, get_furtest_surface_point, check_if_closed_surface
 from vampy.automatedPreprocessing.simulate import run_simulation
 from vampy.automatedPreprocessing.visualize import visualize_model
 
 from fsipy.automatedPreprocessing.preprocessing_common import generate_mesh, distance_to_spheres_solid_thickness, \
-    dist_sphere_spheres, convert_xml_mesh_to_hdf5
+    dist_sphere_spheres, convert_xml_mesh_to_hdf5, convert_vtu_mesh_to_xdmf
 
 
 def run_pre_processing(input_model, verbose_print, smoothing_method, smoothing_factor, smoothing_iterations,
@@ -89,6 +89,7 @@ def run_pre_processing(input_model, verbose_print, smoothing_method, smoothing_f
     file_name_surface_name = base_path + "_remeshed_surface.vtp"
     file_name_xml_mesh = base_path + ".xml"
     file_name_vtu_mesh = base_path + ".vtu"
+    file_name_xdmf_mesh = base_path + ".xdmf"
     region_centerlines = None
 
     if remove_all:
@@ -100,7 +101,7 @@ def run_pre_processing(input_model, verbose_print, smoothing_method, smoothing_f
             file_name_probe_points,
             file_name_voronoi, file_name_voronoi_smooth, file_name_voronoi_surface, file_name_surface_smooth,
             file_name_model_flow_ext, file_name_clipped_model, file_name_flow_centerlines, file_name_surface_name,
-            file_name_xml_mesh, file_name_vtu_mesh,
+            file_name_xml_mesh, file_name_vtu_mesh, file_name_xdmf_mesh,
         ]
         for file in files_to_remove:
             if path.exists(file):
@@ -424,17 +425,14 @@ def run_pre_processing(input_model, verbose_print, smoothing_method, smoothing_f
                                                    solid_thickness,
                                                    solid_thickness_parameters)
 
-        write_mesh(compress_mesh, file_name_surface_name, file_name_vtu_mesh, file_name_xml_mesh,
-                   mesh, remeshed_surface)
-
-        # Add .gz to XML mesh file if compressed
-        if compress_mesh:
-            file_name_xml_mesh = file_name_xml_mesh + ".gz"
+        # Write mesh in VTU format
+        write_polydata(remeshed_surface, file_name_surface_name)
+        write_polydata(mesh, file_name_vtu_mesh)
     else:
         mesh = read_polydata(file_name_vtu_mesh)
 
-    print("--- Converting XML mesh to HDF5\n")
-    convert_xml_mesh_to_hdf5(file_name_xml_mesh)
+    # Convert VTU mesh to XDMF
+    convert_vtu_mesh_to_xdmf(file_name_vtu_mesh, file_name_xdmf_mesh)
 
     network, probe_points = setup_model_network(centerlines, file_name_probe_points, region_center, verbose_print)
 

src/fsipy/automatedPreprocessing/preprocessing_common.py

@@ -4,6 +4,7 @@
 from vtk import vtkPolyData
 from dolfin import Mesh, MeshFunction, File, HDF5File
 from pathlib import Path
+import meshio
 
 
 # Global array names
@@ -203,3 +204,47 @@ def convert_xml_mesh_to_hdf5(file_name_xml_mesh: str, scaling_factor: float = 0.
     hdf.write(mesh, "/mesh")
     hdf.write(boundaries, "/boundaries")
     hdf.write(domains, "/domains")
+
+
+def convert_vtu_mesh_to_xdmf(file_name_vtu_mesh: str, file_name_xdmf_mesh: str) -> None:
+    """
+    Convert a VTU mesh to XDMF format using meshio.
+
+    Args:
+        file_name_vtu_mesh (str): Path to the input VTU mesh file.
+        file_name_xdmf_mesh (str): Path to the output XDMF file.
+    """
+    print("--- Converting VTU mesh to XDMF")
+
+    # Load the VTU mesh
+    vtu_mesh = meshio.read(file_name_vtu_mesh)
+
+    # Extract cell data
+    tetra_data = vtu_mesh.cell_data_dict.get("CellEntityIds", {}).get("tetra", None)
+    triangle_data = vtu_mesh.cell_data_dict.get("CellEntityIds", {}).get("triangle", None)
+
+    # Extract cell types and data
+    tetra_cells = None
+    triangle_cells = None
+    for cell in vtu_mesh.cells:
+        if cell.type == "tetra":
+            tetra_cells = cell.data
+        elif cell.type == "triangle":
+            triangle_cells = cell.data
+
+    # Create mesh objects
+    tetra_mesh = meshio.Mesh(points=vtu_mesh.points, cells={"tetra": tetra_cells},
+                             cell_data={"CellEntityIds": [tetra_data]})
+    triangle_mesh = meshio.Mesh(points=vtu_mesh.points, cells=[("triangle", triangle_cells)],
+                                cell_data={"CellEntityIds": [triangle_data]})
+
+    # Define Path objects
+    tetra_xdmf_path = Path(file_name_xdmf_mesh)
+    triangle_xdmf_path = tetra_xdmf_path.with_name(tetra_xdmf_path.stem + '_triangle.xdmf')
+
+    # Write the VTU mesh to XDMF format
+    meshio.write(tetra_xdmf_path, tetra_mesh)
+    meshio.write(triangle_xdmf_path, triangle_mesh)
+
+    print(f"Tetra mesh XDMF file written to: {tetra_xdmf_path}")
+    print(f"Triangle mesh XDMF file written to: {triangle_xdmf_path}\n")