From 849d66400b6dff53ad44cdb00be7789ba3a5d198 Mon Sep 17 00:00:00 2001 From: Henrik Finsberg Date: Tue, 19 Aug 2025 11:01:55 +0200 Subject: [PATCH 1/3] Change logic for converting mesh from wildmeshing to dolfinx - extract facets on the interface rather than using the data from tracked surfaces coming from wildmeshinh --- src/mri2mesh/mesh/__init__.py | 11 ++- src/mri2mesh/mesh/basic.py | 150 ++++++++++++++++++++++++---------- 2 files changed, 112 insertions(+), 49 deletions(-) diff --git a/src/mri2mesh/mesh/__init__.py b/src/mri2mesh/mesh/__init__.py index 9650ae2..c24bbbe 100644 --- a/src/mri2mesh/mesh/__init__.py +++ b/src/mri2mesh/mesh/__init__.py @@ -27,6 +27,10 @@ def add_mesh_parser(parser: argparse.ArgumentParser) -> None: convert_parser = subparsers.add_parser("convert", help="Convert mesh to dolfinx") convert_parser.add_argument("mesh_dir", type=Path, help="Directory containing mesh files") + convert_parser.add_argument( + "--extract-facet-tags", action="store_true", help="Extract facet tags" + ) + convert_parser.add_argument("--extract-submesh", action="store_true", help="Extract submesh") idealized_parser = subparsers.add_parser( "idealized", @@ -41,12 +45,7 @@ def dispatch(command, args: dict[str, typing.Any]) -> int: basic.generate_sameple_config(**args) elif command == "create": - mesh_dir = basic.create_mesh_from_config(**args) - try: - basic.convert_mesh_dolfinx(mesh_dir=mesh_dir) - except ImportError: - logger.debug("dolfinx not installed, skipping conversion to dolfinx") - + basic.create_mesh_from_config(**args) elif command == "convert": basic.convert_mesh_dolfinx(**args) diff --git a/src/mri2mesh/mesh/basic.py b/src/mri2mesh/mesh/basic.py index d996e82..59d076b 100644 --- a/src/mri2mesh/mesh/basic.py +++ b/src/mri2mesh/mesh/basic.py @@ -156,73 +156,102 @@ def create_mesh( manifold_surface=manifold_surface, correct_surface_orientation=True, ) - coords, connections = tetra.get_tracked_surfaces() tetra_mesh = meshio.Mesh( point_array, [("tetra", cell_array)], cell_data={"cell_tags": [marker.ravel()]} ) - tetra_mesh_pv = pv.from_meshio(tetra_mesh).clean() + + tetra_mesh_pv = pv.from_meshio(tetra_mesh) pv.save_meshio(outdir / "tetra_mesh.xdmf", tetra_mesh_pv) - for i, coord in enumerate(coords): - np.save(outdir / f"coords_{i}.npy", coord) + np.save(outdir / "point_array.npy", point_array) + np.save(outdir / "cell_array.npy", cell_array) + np.save(outdir / "marker.npy", marker) + + # coords, connections = tetra.get_tracked_surfaces() + # for i, coord in enumerate(coords): + # np.save(outdir / f"coords_{i}.npy", coord) + + # for i, conn in enumerate(connections): + # np.save(outdir / f"connections_{i}.npy", conn) - for i, conn in enumerate(connections): - np.save(outdir / f"connections_{i}.npy", conn) + # cell_tags.find + # Compute incident with facets + # Intersection exterior facets -def convert_mesh_dolfinx(mesh_dir: Path): +def convert_mesh_dolfinx( + mesh_dir: Path, extract_facet_tags: bool = False, extract_submesh: bool = False +): logger.info("Converting mesh to dolfinx in %s", mesh_dir) from mpi4py import MPI - from scipy.spatial.distance import cdist import dolfinx + import basix + import ufl - threshold = 1.0 - fdim = 2 - - coords = [] - for path in sorted(mesh_dir.glob("coords_*.npy"), key=lambda x: int(x.stem.split("_")[-1])): - coords.append(np.load(path)) - logger.debug(f"Found {len(coords)} coordinates") + point_array = np.load(mesh_dir / "point_array.npy") + cell_array = np.load(mesh_dir / "cell_array.npy") + marker = np.load(mesh_dir / "marker.npy") - connections = [] - for path in sorted( - mesh_dir.glob("connections_*.npy"), key=lambda x: int(x.stem.split("_")[-1]) - ): - connections.append(np.load(path)) - logger.debug(f"Found {len(connections)} connections") + comm = MPI.COMM_WORLD + mesh = dolfinx.mesh.create_mesh( + comm, + cell_array.astype(np.int64), + point_array, + ufl.Mesh(basix.ufl.element("Lagrange", "tetrahedron", 1, shape=(3,))), + ) + tdim = mesh.topology.dim + fdim = tdim - 1 + local_entities, local_values = dolfinx.io.gmshio.distribute_entity_data( + mesh, + tdim, + cell_array.astype(np.int64), + marker.flatten().astype(np.int32), + ) + adj = dolfinx.graph.adjacencylist(local_entities) + cell_tags = dolfinx.mesh.meshtags_from_entities( + mesh, + tdim, + adj, + local_values.astype(np.int32, copy=False), + ) + if not extract_facet_tags: + logger.debug("Save files") + with dolfinx.io.XDMFFile(comm, mesh_dir / "mesh_full.xdmf", "w") as xdmf: + xdmf.write_mesh(mesh) + xdmf.write_meshtags(cell_tags, mesh.geometry) - assert len(connections) == len(coords) + return - logger.debug("Loading mesh") - comm = MPI.COMM_WORLD - with dolfinx.io.XDMFFile(comm, mesh_dir / "tetra_mesh.xdmf", "r") as xdmf: - mesh = xdmf.read_mesh(name="Grid") - cell_tags = xdmf.read_meshtags(mesh, name="Grid") + mesh.topology.create_connectivity(tdim - 1, tdim) - logger.debug("Mesh loaded") + # FIXME: Here we just add hard coded values for now. This should be fixed in the future - facets = [] + entities = [] values = [] - for i, coord in enumerate(coords, start=1): - logger.debug(f"Processing coord {i}") + # 1 = Parenchyma + PARENCHYMA = 1 - def locator(x): - # Find the distance to all coordinates - distances = cdist(x.T, coord) - # And return True is they are close - return np.any(distances < threshold, axis=1) + cells = cell_tags.find(PARENCHYMA) + exterior_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology) + incident_facets = dolfinx.mesh.compute_incident_entities(mesh.topology, cells, tdim, tdim - 1) + exterior_facets_marker = np.intersect1d(incident_facets, exterior_facets) + values.append(np.full(exterior_facets_marker.shape[0], PARENCHYMA, dtype=np.int32)) + entities.append(exterior_facets_marker) - f = dolfinx.mesh.locate_entities_boundary(mesh, dim=fdim, marker=locator) - v = np.full(f.shape[0], i, dtype=np.int32) - facets.append(f) - values.append(v) + VENTRICLES = 3 + all_cell_tags = np.unique(cell_tags.values) + cell_not_ventricles = np.setdiff1d(all_cell_tags, [VENTRICLES]) + import scifem + + interface_entities = scifem.mesh.find_interface(cell_tags, [VENTRICLES], cell_not_ventricles) + entities.append(interface_entities) + values.append(np.full(interface_entities.shape[0], VENTRICLES, dtype=np.int32)) - logger.debug("Create meshtags") facet_tags = dolfinx.mesh.meshtags( mesh, fdim, - np.hstack(facets), + np.hstack(entities), np.hstack(values), ) facet_tags.name = "facet_tags" @@ -230,9 +259,44 @@ def locator(x): mesh.name = "mesh" logger.debug("Save files") - with dolfinx.io.XDMFFile(comm, mesh_dir / "mesh.xdmf", "w") as xdmf: + with dolfinx.io.XDMFFile(comm, mesh_dir / "mesh_full.xdmf", "w") as xdmf: xdmf.write_mesh(mesh) xdmf.write_meshtags(facet_tags, mesh.geometry) xdmf.write_meshtags(cell_tags, mesh.geometry) logger.info("Mesh saved to %s", mesh_dir / "mesh.xdmf") + + if not extract_submesh: + return + submesh_data = scifem.mesh.extract_submesh(mesh, cell_tags, cell_not_ventricles) + + # Transfer the facet tags to the submesh + submesh_data.domain.topology.create_connectivity(2, 3) + facet_tags_submesh, sub_to_parent_entity_map = scifem.mesh.transfer_meshtags_to_submesh( + facet_tags, + # geo.facet_tags, # If available + submesh_data.domain, + vertex_to_parent=submesh_data.vertex_map, + cell_to_parent=submesh_data.cell_map, + ) + + np.save(mesh_dir / "sub_to_parent_entity_map.npy", sub_to_parent_entity_map) + np.save(mesh_dir / "vertex_map.npy", submesh_data.vertex_map) + np.save(mesh_dir / "cell_map.npy", submesh_data.cell_map) + + # Remove overflow values + keep_indices = facet_tags_submesh.values > 0 + facet_tags_new = dolfinx.mesh.meshtags( + submesh_data.domain, + fdim, + facet_tags_submesh.indices[keep_indices], + facet_tags_submesh.values[keep_indices], + ) + + facet_tags_new.name = "facet_tags" + submesh_data.cell_tag.name = "cell_tags" + + with dolfinx.io.XDMFFile(comm, mesh_dir / "mesh.xdmf", "w") as xdmf: + xdmf.write_mesh(submesh_data.domain) + xdmf.write_meshtags(submesh_data.cell_tag, submesh_data.domain.geometry) + xdmf.write_meshtags(facet_tags_new, submesh_data.domain.geometry) From de93136ac6b1e9c793551266235945741c08b553 Mon Sep 17 00:00:00 2001 From: Henrik Finsberg Date: Tue, 19 Aug 2025 11:55:48 +0200 Subject: [PATCH 2/3] Attempt to fix tests --- examples/cli.ipynb | 36 ++++++++++++++++++++++++++++++++---- src/mri2mesh/mesh/basic.py | 9 +++++---- tests/test_cli.py | 13 ++++--------- 3 files changed, 41 insertions(+), 17 deletions(-) diff --git a/examples/cli.ipynb b/examples/cli.ipynb index 29222d9..4c03a62 100644 --- a/examples/cli.ipynb +++ b/examples/cli.ipynb @@ -283,7 +283,7 @@ "id": "27", "metadata": {}, "source": [ - "Let us have a look at the mesh" + "Now lets convert the mesh to dolfinx" ] }, { @@ -292,6 +292,34 @@ "id": "28", "metadata": {}, "outputs": [], + "source": [ + "!mri2mesh mesh convert idealized-brain" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "29", + "metadata": {}, + "outputs": [], + "source": [ + "!ls idealized-brain/" + ] + }, + { + "cell_type": "markdown", + "id": "30", + "metadata": {}, + "source": [ + "Let us have a look at the mesh" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "31", + "metadata": {}, + "outputs": [], "source": [ "from mpi4py import MPI\n", "import dolfinx\n", @@ -304,7 +332,7 @@ { "cell_type": "code", "execution_count": null, - "id": "29", + "id": "32", "metadata": {}, "outputs": [], "source": [ @@ -326,7 +354,7 @@ { "cell_type": "code", "execution_count": null, - "id": "30", + "id": "33", "metadata": {}, "outputs": [], "source": [] @@ -348,7 +376,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.12.9" + "version": "3.12.11" } }, "nbformat": 4, diff --git a/src/mri2mesh/mesh/basic.py b/src/mri2mesh/mesh/basic.py index 59d076b..6334e37 100644 --- a/src/mri2mesh/mesh/basic.py +++ b/src/mri2mesh/mesh/basic.py @@ -217,7 +217,7 @@ def convert_mesh_dolfinx( ) if not extract_facet_tags: logger.debug("Save files") - with dolfinx.io.XDMFFile(comm, mesh_dir / "mesh_full.xdmf", "w") as xdmf: + with dolfinx.io.XDMFFile(comm, mesh_dir / "mesh.xdmf", "w") as xdmf: xdmf.write_mesh(mesh) xdmf.write_meshtags(cell_tags, mesh.geometry) @@ -225,7 +225,7 @@ def convert_mesh_dolfinx( mesh.topology.create_connectivity(tdim - 1, tdim) - # FIXME: Here we just add hard coded values for now. This should be fixed in the future + # FIXME: Here we just add hard coded values for now. This should be fixed in the future. entities = [] values = [] @@ -259,12 +259,13 @@ def convert_mesh_dolfinx( mesh.name = "mesh" logger.debug("Save files") - with dolfinx.io.XDMFFile(comm, mesh_dir / "mesh_full.xdmf", "w") as xdmf: + meshname = "mesh_full.xdmf" if extract_submesh else "mesh.xdmf" + with dolfinx.io.XDMFFile(comm, mesh_dir / meshname, "w") as xdmf: xdmf.write_mesh(mesh) xdmf.write_meshtags(facet_tags, mesh.geometry) xdmf.write_meshtags(cell_tags, mesh.geometry) - logger.info("Mesh saved to %s", mesh_dir / "mesh.xdmf") + logger.info("Mesh saved to %s", mesh_dir / meshname) if not extract_submesh: return diff --git a/tests/test_cli.py b/tests/test_cli.py index 054a41f..2d2ceef 100644 --- a/tests/test_cli.py +++ b/tests/test_cli.py @@ -14,14 +14,9 @@ def test_labels_cli(name, capsys): def test_mesh_idealized_cli(tmp_path): cli.main(["mesh", "idealized", "-o", str(tmp_path)]) for name in [ - "connections_3.npy", - "coords_0.npy", - "coords_1.npy", - "connections_2.npy", - "connections_0.npy", - "coords_3.npy", - "coords_2.npy", - "connections_1.npy", + "point_array.npy", + "cell_array.npy", + "marker.npy", "mesh.h5", "ventricles.ply", "skull.ply", @@ -34,4 +29,4 @@ def test_mesh_idealized_cli(tmp_path): "V34.ply", "mesh_params.json", ]: - assert (tmp_path / name).exists() + assert (tmp_path / name).exists(), name From c12f04ad66c43209f7dd3e09386950e8031724ef Mon Sep 17 00:00:00 2001 From: Henrik Finsberg Date: Tue, 19 Aug 2025 12:14:06 +0200 Subject: [PATCH 3/3] Add name to cell tags before saving --- src/mri2mesh/mesh/basic.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/mri2mesh/mesh/basic.py b/src/mri2mesh/mesh/basic.py index 6334e37..f0964de 100644 --- a/src/mri2mesh/mesh/basic.py +++ b/src/mri2mesh/mesh/basic.py @@ -215,6 +215,7 @@ def convert_mesh_dolfinx( adj, local_values.astype(np.int32, copy=False), ) + cell_tags.name = "cell_tags" if not extract_facet_tags: logger.debug("Save files") with dolfinx.io.XDMFFile(comm, mesh_dir / "mesh.xdmf", "w") as xdmf: @@ -255,7 +256,7 @@ def convert_mesh_dolfinx( np.hstack(values), ) facet_tags.name = "facet_tags" - cell_tags.name = "cell_tags" + mesh.name = "mesh" logger.debug("Save files")