FIX: fix volume meshing of cavities

src/ansys/heart/core/helpers/fluenthdf5.py

@@ -426,7 +426,9 @@ def _convert_interior_faces_to_tetrahedrons2(self) -> Tuple[np.ndarray, np.ndarr
         return tetrahedrons, cell_ids
 
     # NOTE: no typehint due to lazy import of pyvista
-    def _to_vtk(self, add_cells: bool = True, add_faces: bool = False):
+    def _to_vtk(
+        self, add_cells: bool = True, add_faces: bool = False, remove_interior_faces: bool = False
+    ):
         """Convert mesh to vtk unstructured grid or polydata.
 
         Parameters
@@ -471,11 +473,16 @@ def _to_vtk(self, add_cells: bool = True, add_faces: bool = False):
 
         if add_faces:
             # add faces.
+            face_zones = self.face_zones
+
+            if remove_interior_faces:
+                face_zones = [fz for fz in face_zones if "interior" not in fz.name]
+
             grid_faces = pv.UnstructuredGrid()
             grid_faces.nodes = self.nodes
 
-            face_zone_ids = np.concatenate([[fz.id] * fz.faces.shape[0] for fz in self.face_zones])
-            faces = np.array(np.concatenate([fz.faces for fz in self.face_zones]), dtype=int)
+            face_zone_ids = np.concatenate([[fz.id] * fz.faces.shape[0] for fz in face_zones])
+            faces = np.array(np.concatenate([fz.faces for fz in face_zones]), dtype=int)
             faces = np.hstack([np.ones((faces.shape[0], 1), dtype=int) * 3, faces])
 
             grid_faces = pv.UnstructuredGrid(

src/ansys/heart/core/models.py

@@ -616,14 +616,8 @@ def mesh_volume(
 
         return
 
-    def _mesh_fluid_volume(self, remesh_caps: bool = True):
-        """Generate a volume mesh of the cavities.
-
-        Parameters
-        ----------
-        remesh_caps : bool, optional
-            Flag indicating whether to remesh the caps of each cavity, by default True
-        """
+    def _mesh_fluid_volume(self):
+        """Generate a volume mesh of the cavities."""
         # get all relevant boundaries for the fluid cavities:
         substrings_include = ["endocardium", "valve-plane", "septum"]
         substrings_include_re = "|".join(substrings_include)
@@ -651,38 +645,45 @@ def _mesh_fluid_volume(self, remesh_caps: bool = True):
 
         LOGGER.info("Meshing fluid cavities...")
 
+        # get list of fluid cavities
         # mesh the fluid cavities
-        fluid_mesh = mesher.mesh_fluid_cavities(
-            boundaries_fluid, caps, self.workdir, remesh_caps=remesh_caps
-        )
+        cavity_surfaces = [
+            self.mesh.get_surface(part.cavity.surface.id) for part in self.parts if part.cavity
+        ]
+        # remove caps.
+        cavity_surfaces = [
+            SurfaceMesh(cs.threshold((0, 0), "_cap_id").extract_surface(), name=cs.name)
+            for cs in cavity_surfaces
+        ]
+        fluid_mesh = mesher._mesh_fluid_cavities(cavity_surfaces, self.workdir, mesh_size=1)
 
-        LOGGER.info(f"Meshed {len(fluid_mesh.cell_zones)} fluid regions...")
+        # LOGGER.info(f"Meshed {len(fluid_mesh.cell_zones)} fluid regions...")
 
         # add part-ids
-        cz_ids = np.sort([cz.id for cz in fluid_mesh.cell_zones])
+        # cz_ids = np.sort([cz.id for cz in fluid_mesh.cell_zones])
 
         # TODO: this offset is arbitrary.
-        offset = 10000
-        new_ids = np.arange(cz_ids.shape[0]) + offset
-        czid_to_pid = {cz_id: new_ids[ii] for ii, cz_id in enumerate(cz_ids)}
+        # offset = 10000
+        # new_ids = np.arange(cz_ids.shape[0]) + offset
+        # czid_to_pid = {cz_id: new_ids[ii] for ii, cz_id in enumerate(cz_ids)}
 
-        for cz in fluid_mesh.cell_zones:
-            cz.id = czid_to_pid[cz.id]
+        # for cz in fluid_mesh.cell_zones:
+        #     cz.id = czid_to_pid[cz.id]
 
-        fluid_mesh._fix_negative_cells()
-        fluid_mesh_vtk = fluid_mesh._to_vtk(add_cells=True, add_faces=False)
+        # fluid_mesh._fix_negative_cells()
+        # fluid_mesh_vtk = fluid_mesh._to_vtk(add_cells=True, add_faces=False)
 
-        fluid_mesh_vtk.cell_data["part-id"] = fluid_mesh_vtk.cell_data["cell-zone-ids"]
+        # fluid_mesh_vtk.cell_data["part-id"] = fluid_mesh_vtk.cell_data["cell-zone-ids"]
 
-        boundaries = [
-            SurfaceMesh(name=fz.name, triangles=fz.faces, nodes=fluid_mesh.nodes, id=fz.id)
-            for fz in fluid_mesh.face_zones
-            if "interior" not in fz.name
-        ]
+        # boundaries = [
+        #     SurfaceMesh(name=fz.name, triangles=fz.faces, nodes=fluid_mesh.nodes, id=fz.id)
+        #     for fz in fluid_mesh.face_zones
+        #     if "interior" not in fz.name
+        # ]
 
-        self.fluid_mesh = Mesh(fluid_mesh_vtk)
-        for boundary in boundaries:
-            self.fluid_mesh.add_surface(boundary, boundary.id, boundary.name)
+        self.fluid_mesh = Mesh(fluid_mesh)
+        # for boundary in boundaries:
+        #     self.fluid_mesh.add_surface(boundary, boundary.id, boundary.name)
 
         return
 
@@ -1412,6 +1413,7 @@ def _assign_cavities_to_parts(self) -> None:
 
                 # create cap: NOTE, mostly for compatibility. Could simplify further
                 cap_mesh = SurfaceMesh(patch.clean(), name=cap_name, id=ii + idoffset)
+                cap_mesh.cell_data["_cap_id"] = float(ii + idoffset)
 
                 # Add cap to main mesh.
                 self.mesh.add_surface(cap_mesh, id=cap_mesh.id, name=cap_name)
@@ -1431,7 +1433,7 @@ def _assign_cavities_to_parts(self) -> None:
             # ! and we can only use this for getting some meta-data such as volume
             # ! also it is not updated dynamically.
             # merge patches into cavity surface.
-            surface.cell_data["_cap_id"] = 0
+            surface.cell_data["_cap_id"] = float(0)
             surface_cavity = SurfaceMesh(pv.merge([surface] + [cap._mesh for cap in part.caps]))
             surface_cavity.name = surface.name
             surface_cavity.id = surface.id

src/ansys/heart/preprocessor/mesher.py

@@ -232,101 +232,6 @@ def _wrap_part(session: MeshingSession, boundary_names: list, wrapped_part_name:
     return wrapped_face_zone_names
 
 
-def mesh_fluid_cavities(
-    fluid_boundaries: List[SurfaceMesh],
-    caps: List[SurfaceMesh],
-    workdir: str,
-    remesh_caps: bool = True,
-) -> FluentMesh:
-    """Mesh the fluid cavities.
-
-    Parameters
-    ----------
-    fluid_boundaries : List[SurfaceMesh]
-        List of fluid boundaries used for meshing.
-    caps : List[SurfaceMesh]
-        List of caps that close each of the cavities.
-    workdir : str
-        Working directory
-    remesh_caps : bool, optional
-        Flag indicating whether to remesh the caps, by default True
-
-    Returns
-    -------
-    Path
-        Path to the .msh.h5 volume mesh.
-    """
-    if _uses_container:
-        mounted_volume = pyfluent.EXAMPLES_PATH
-        work_dir_meshing = os.path.join(mounted_volume, "tmp_meshing-fluid")
-    else:
-        work_dir_meshing = os.path.join(workdir, "meshing-fluid")
-
-    if not os.path.isdir(work_dir_meshing):
-        os.makedirs(work_dir_meshing)
-    else:
-        files = glob.glob(os.path.join(work_dir_meshing, "*.stl"))
-        for f in files:
-            os.remove(f)
-
-    # write all boundaries
-    for b in fluid_boundaries:
-        filename = os.path.join(work_dir_meshing, b.name.lower() + ".stl")
-        b.save(filename)
-        add_solid_name_to_stl(filename, b.name.lower(), file_type="binary")
-
-    for c in caps:
-        filename = os.path.join(work_dir_meshing, c.name.lower() + ".stl")
-        c.save(filename)
-        add_solid_name_to_stl(filename, c.name.lower(), file_type="binary")
-
-    session = _get_fluent_meshing_session(work_dir_meshing)
-
-    # import all stls
-    if _uses_container:
-        # NOTE: when using a Fluent container visible files
-        # will be in /mnt/pyfluent. So need to use relative paths
-        # or replace dirname by /mnt/pyfluent as prefix
-        work_dir_meshing = "/mnt/pyfluent/meshing"
-    session.tui.file.import_.cad(f"no {work_dir_meshing} *.stl")
-
-    # merge objects
-    session.tui.objects.merge("'(*)", "model-fluid")
-
-    # fix duplicate nodes
-    session.tui.diagnostics.face_connectivity.fix_free_faces("objects '(*)")
-
-    # set size field
-    session.tui.size_functions.set_global_controls(1, 1, 1.2)
-    session.tui.scoped_sizing.compute("yes")
-
-    # remesh all caps
-    if remesh_caps:
-        session.tui.boundary.remesh.remesh_constant_size("(cap_*)", "()", 40, 20, 1, "yes")
-
-    # convert to mesh object
-    session.tui.objects.change_object_type("(*)", "mesh", "yes")
-
-    # compute volumetric regions
-    session.tui.objects.volumetric_regions.compute("model-fluid")
-
-    # mesh volume
-    session.tui.mesh.auto_mesh("model-fluid")
-
-    # clean up
-    session.tui.objects.delete_all_geom()
-    session.tui.objects.delete_unreferenced_faces_and_edges()
-
-    # write
-    file_path_mesh = os.path.join(workdir, "fluid-mesh.msh.h5")
-    session.tui.file.write_mesh(file_path_mesh)
-
-    mesh = FluentMesh(file_path_mesh)
-    mesh.load_mesh()
-
-    return mesh
-
-
 def mesh_from_manifold_input_model(
     model: _InputModel,
     workdir: Union[str, Path],
@@ -765,5 +670,125 @@ def mesh_from_non_manifold_input_model(
     return new_mesh
 
 
+def _mesh_fluid_cavities(
+    cavity_boundaries: list[SurfaceMesh], workdir: str, mesh_size: float = 1.0
+) -> pv.UnstructuredGrid:
+    """Mesh the caps of each fluid cavity with uniformly.
+
+    Parameters
+    ----------
+    cavity_boundaries : List[SurfaceMesh]
+        List of cavity boundaries used for meshing.
+    workdir : str
+        Working directory
+    mesh_size : float
+        Mesh size
+
+    Returns
+    -------
+    pv.UnstructuredGrid
+        Unstructured grid with fluid mesh.
+    """
+    # Use the following tgrid api utility for each cavity:
+    # Consequently need to associate each "patch" with a cap. E.g. by centroid?
+    #
+    # (tgapi-util-fill-holes-in-face-zone-list '(face-zone-list) max-hole-edges)
+    if _uses_container:
+        mounted_volume = pyfluent.EXAMPLES_PATH
+        work_dir_meshing = os.path.join(mounted_volume, "tmp_meshing-fluid")
+    else:
+        work_dir_meshing = os.path.join(workdir, "meshing-fluid")
+
+    if not os.path.isdir(work_dir_meshing):
+        os.makedirs(work_dir_meshing)
+    else:
+        files = glob.glob(os.path.join(work_dir_meshing, "*.stl"))
+        for f in files:
+            os.remove(f)
+
+    # write all boundaries
+    for b in cavity_boundaries:
+        filename = os.path.join(work_dir_meshing, b.name.lower() + ".stl")
+        b.save(filename)
+        add_solid_name_to_stl(filename, b.name.lower(), file_type="binary")
+
+    session = _get_fluent_meshing_session(work_dir_meshing)
+
+    # import all stls
+    if _uses_container:
+        # NOTE: when using a Fluent container visible files
+        # will be in /mnt/pyfluent. So need to use relative paths
+        # or replace dirname by /mnt/pyfluent as prefix
+        work_dir_meshing = "/mnt/pyfluent/meshing"
+
+    session.tui.file.import_.cad(f"no {work_dir_meshing} *.stl")
+
+    # set size field
+    session.tui.size_functions.set_global_controls(mesh_size, mesh_size, 1.2)
+    session.tui.scoped_sizing.compute("yes")
+
+    # create caps
+    # (tgapi-util-fill-holes-in-face-zone-list '(face-zone-list) max-hole-edges)
+    # convert all to mesh object
+    session.tui.objects.change_object_type("'(*)", "mesh", "yes")
+    for cavity_boundary in cavity_boundaries:
+        cavity_name = "-".join(cavity_boundary.name.split()).lower()
+        session.scheme_eval.scheme_eval(
+            f"(tgapi-util-fill-holes-in-face-zone-list '({cavity_name}) 1000)"
+        )
+        # merge unreferenced with cavity
+        # (get-unreferenced-face-zones)
+        patch_ids = session.scheme_eval.scheme_eval("(get-unreferenced-face-zones)")
+        session.tui.objects.create(
+            f"{cavity_name}-patches",
+            "fluid",
+            3,
+            "({0})".format(" ".join([str(patch_id) for patch_id in patch_ids])),
+            "()",
+            "mesh",
+            "yes",
+        )
+        # merge objects
+        session.tui.objects.merge(f"'({cavity_name}*)")
+
+    # find overlapping pairs
+    # (tgapi-util-get-overlapping-face-zones "*patch*" 0.1 0.1 )
+    # try:
+    #     overlapping_pairs = session.scheme_eval.scheme_eval(
+    #         '(tgapi-util-get-overlapping-face-zones "*patch*" 0.1 0.1 )'
+    #     )
+    #     for pair in overlapping_pairs:
+    #         # remove first from pair
+    #         session.tui.boundary.manage.delete(f"'({pair[0]})", "yes")
+    # except:
+    #     LOGGER.debug("No overlapping face zones.")
+
+    # merge mesh objects
+    mesh_object_names = session.scheme_eval.scheme_eval("(get-objects-of-type 'mesh)")
+    if len(mesh_object_names) > 1:
+        session.tui.objects.merge("'(*)", "fluid-zones")
+    else:
+        session.tui.objects.rename_object(mesh_object_names[0].str, "fluid-zones")
+
+    # merge nodes
+    session.tui.boundary.merge_nodes("'(*)", "'(*)", "no", "no , ,")
+    # mesh volumes
+    session.tui.objects.volumetric_regions.compute("fluid-zones", "no")
+    session.tui.mesh.auto_mesh("fluid-zones", "yes", "pyr", "tet", "yes")
+
+    session.tui.objects.delete_all_geom()
+
+    file_path_mesh = os.path.join(workdir, "fluid-mesh.msh.h5")
+    session.tui.file.write_mesh(file_path_mesh)
+
+    session.exit()
+
+    # # write
+    mesh = FluentMesh(file_path_mesh)
+    mesh.load_mesh(reconstruct_tetrahedrons=True)
+
+    return mesh._to_vtk(add_cells=True, add_faces=True, remove_interior_faces=True)
+
+
 if __name__ == "__main__":
     LOGGER.info("Protected")