fix: update default parameters for mesh_fix

Set verbose and joincomp parameters to True by default in the mesh_fix function to improve clarity and ensure consistent behavior.

src/mkit/io/_meshio.py

@@ -30,6 +30,12 @@ def load_meshio(filename: StrPath) -> meshio.Mesh:
 def as_meshio(mesh: AnyMesh, **kwargs: Unpack[Attrs]) -> meshio.Mesh:
     match mesh:
         case meshio.Mesh():
+            if point_data := kwargs.get("point_data"):
+                mesh.point_data.update(point_data)
+            if cell_data := kwargs.get("cell_data"):
+                mesh.cell_data.update(cell_data)
+            if field_data := kwargs.get("field_data"):
+                mesh.field_data.update(field_data)
             return mesh
         case pytorch3d.structures.Meshes():
             return pytorch3d_to_meshio(mesh, **kwargs)

src/mkit/ops/mesh_fix.py

@@ -8,8 +8,8 @@
 def mesh_fix(
     mesh: trimesh.Trimesh,
     *,
-    verbose: bool = False,
-    joincomp: bool = False,
+    verbose: bool = True,
+    joincomp: bool = True,
     remove_smallest_components: bool = True,
 ) -> trimesh.Trimesh:
     mesh_pv: pv.PolyData = mkit.io.as_pyvista(mesh)

src/mkit/ops/tetgen.py

@@ -12,10 +12,11 @@
 from mkit.typing import StrPath
 
 
-def tetgen(mesh: meshio.Mesh) -> meshio.Mesh:
+def tetgen(mesh: meshio.Mesh, *, verbose: bool = True) -> meshio.Mesh:
     """
     Args:
         mesh: input mesh
+        holes: (N, 3) float
 
     Returns:
         tetrahedral mesh
@@ -24,10 +25,10 @@ def tetgen(mesh: meshio.Mesh) -> meshio.Mesh:
         tmpdir = pathlib.Path(tmpdir)
         input_file: pathlib.Path = tmpdir / "mesh.smesh"
         save_smesh(input_file, mesh)
-        subprocess.run(
-            ["tetgen", "-p", "-q", "-O", "-z", "-k", "-C", "-V", input_file],
-            check=True,
-        )
+        args: list[str] = ["tetgen", "-p", "-q", "-O", "-z", "-k", "-C"]
+        if verbose:
+            args.append("-V")
+        subprocess.run([*args, input_file], check=True)
         tetra_mesh: meshio.Mesh = meshio.read(tmpdir / "mesh.1.vtk")
         points: npt.NDArray[np.floating] = tetra_mesh.points
         tetra: npt.NDArray[np.intp] = tetra_mesh.get_cells_type("tetra")
@@ -109,7 +110,7 @@ def save_smesh(file: StrPath, mesh: meshio.Mesh) -> None:
         fprint()
         fprint("# Part 3 - hole list")
         fprint("# <# of holes>")
-        holes: npt.NDArray[np.float64] = mesh.field_data.get("holes", [])
+        holes: npt.NDArray[np.floating] = np.asarray(mesh.field_data.get("holes", []))
         fprint(len(holes))
         fprint("# <hole #> <x> <y> <z>")
         for hole_id, hole in enumerate(holes):

tasks/register/Makefile

@@ -33,6 +33,7 @@ $(TEMPLATE_DIR)/01-face.ply: $(TEMPLATE_DIR)/00-face.ply $(srcdir)/template/face
 $(TEMPLATE_DIR)/01-skull.vtu $(TEMPLATE_DIR)/01-mandible.vtu $(TEMPLATE_DIR)/01-maxilla.vtu &: $(TEMPLATE_DIR)/00-skull.ply $(srcdir)/template/skull/process.py
 	python "$(lastword $^)" --output="$(TEMPLATE_DIR)/01-skull.vtu" --mandible="$(TEMPLATE_DIR)/01-mandible.vtu" --maxilla="$(TEMPLATE_DIR)/01-maxilla.vtu" "$<"
 
+# Trim template face according to target face
 ground: $(TEMPLATE_DIR)/00-face.ply $(TARGET_DIR)/01-face.ply $(srcdir)/template/face/ground.py
 	python "$(lastword $^)" "$<" "$(word 2,$^)"
 

tasks/tetgen/src/create.py

@@ -1,16 +1,47 @@
 import pathlib
 from typing import Annotated
 
+import meshio
 import mkit.cli
+import mkit.io
+import mkit.ops.ray
+import mkit.ops.tetgen
+import numpy as np
+import scipy.interpolate
+import trimesh
 import typer
+from numpy import typing as npt
 
 
 def main(
-    face_file: Annotated[pathlib.Path, typer.Option("--face", exists=True)],
-    mandible_file: Annotated[pathlib.Path, typer.Option("--mandible", exists=True)],
-    maxilla_file: Annotated[pathlib.Path, typer.Option("--maxilla", exists=True)],
+    *,
+    output_file: Annotated[
+        pathlib.Path, typer.Option("-o", "--output", dir_okay=False, writable=True)
+    ],
 ) -> None:
-    pass
+    face: trimesh.Trimesh = trimesh.creation.icosphere(radius=1.0)
+    pre_skull: trimesh.Trimesh = trimesh.creation.icosphere(radius=0.5)
+    post_skull: trimesh.Trimesh = pre_skull.copy()
+    post_skull.vertices[post_skull.vertices[:, 1] < 0.0] += [0.0, -0.1, 0.0]
+    tri: trimesh.Trimesh = trimesh.util.concatenate([face, pre_skull])
+    tet: meshio.Mesh = mkit.ops.tetgen.tetgen(
+        mkit.io.as_meshio(
+            tri, field_data={"holes": [mkit.ops.ray.find_inner_point(pre_skull)]}
+        )
+    )
+    closest: npt.NDArray[np.floating]
+    distance: npt.NDArray[np.floating]
+    triangle_id: npt.NDArray[np.integer]
+    closest, distance, triangle_id = pre_skull.nearest.on_surface(tet.points)
+    mask_skull: npt.NDArray[np.bool_] = distance < 1e-6
+    disp: npt.NDArray[np.floating] = np.full(tet.points.shape, np.nan)
+    disp[mask_skull] = scipy.interpolate.griddata(
+        pre_skull.vertices,
+        post_skull.vertices - pre_skull.vertices,
+        tet.points[mask_skull],
+        method="nearest",
+    )
+    mkit.io.save(output_file, tet, point_data={"disp": disp})
 
 
 if __name__ == "__main__":