Generating meshes from STL files#
In this demo, we will show how to generate a mesh from STL files based on MRI images. We use the files from MRI2FEM2: Chapter 3.
We start by importing all the packages required for this task.
from mpi4py import MPI
import wildmeshing
import dolfinx
from pathlib import Path
import json
import numpy as np
import pyvista
import ufl
import basix.ufl
import os
Next, we use the already downloaded stl files from this chapter. We will use WildMeshing to combine these. We start by taking the union of the two choroid plexus files (left and right hemisphere). Next, we embed these within the ventricles. The operations are described as a dictionary structure below
folder = Path(os.environ["WILDFENICS_DATA_PATH"])
assert folder.exists(), "Could not find stl files"
tree = {
"operation": "union",
"left": str((folder / "dura.final.stl").absolute().as_posix()),
"right": {
"operation": "union",
"left": str((folder / "white.final.stl").absolute().as_posix()),
"right": {
"operation": "union",
"left": str((folder / "ventricles.final.stl").absolute().as_posix()),
"right": {
"operation": "union",
"left": str(
(folder / "rhchoroid.plexus.final.stl").absolute().as_posix()
),
"right": str(
(folder / "lhchoroid.plexus.final.stl").absolute().as_posix()
),
},
},
},
}
Next, we can set up a 3D meshing instance and load in the instruction-set
tetra = wildmeshing.Tetrahedralizer(stop_quality=100, max_its=30, edge_length_r=0.1)
tetra.load_csg_tree(json.dumps(tree))
TBB threads 4
o-[I/O ] Loading file /root/data/mri2femii-chp2-dataset/Gonzo/output/dura.final.stl...
Error: Could not load file: /root/data/mri2femii-chp2-dataset/Gonzo/output/dura.final.stl
[2025-03-03 13:16:30.235] [float-tetwild] [error] unable to open /root/data/mri2femii-chp2-dataset/Gonzo/output/dura.final.stl file
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[3], line 2
1 tetra = wildmeshing.Tetrahedralizer(stop_quality=100, max_its=30, edge_length_r=0.1)
----> 2 tetra.load_csg_tree(json.dumps(tree))
ValueError: Invalid mesh path in the json
Then, we can generate the mesh
# Create mesh
tetra.tetrahedralize()
[2025-02-17 08:58:13.832] [float-tetwild] [info] remove duplicates:
[2025-02-17 08:58:13.832] [float-tetwild] [info] #v: 141972 -> 23655
[2025-02-17 08:58:13.832] [float-tetwild] [info] #f: 47324 -> 47324
collapsing 0.897446
swapping 0.015215
[2025-02-17 08:58:14.751] [float-tetwild] [info] remove duplicates:
[2025-02-17 08:58:14.751] [float-tetwild] [info] #v: 6583 -> 6583
[2025-02-17 08:58:14.751] [float-tetwild] [info] #f: 13180 -> 13180
[2025-02-17 08:58:14.751] [float-tetwild] [info] #v = 6583
[2025-02-17 08:58:14.751] [float-tetwild] [info] #f = 13179
#boundary_e1 = 1
#boundary_e2 = 2
[2025-02-17 08:58:14.762] [float-tetwild] [info] preprocessing 0.960763s
[2025-02-17 08:58:14.762] [float-tetwild] [info]
[2025-02-17 08:58:14.817] [float-tetwild] [info] #v = 8433
[2025-02-17 08:58:14.817] [float-tetwild] [info] #t = 50157
[2025-02-17 08:58:14.817] [float-tetwild] [info] tetrahedralizing 0.055105s
[2025-02-17 08:58:14.817] [float-tetwild] [info]
[2025-02-17 08:58:14.817] [float-tetwild] [info] triangle insertion start, #f = 13179, #v = 8433, #t = 50157
[2025-02-17 08:58:14.831] [float-tetwild] [info] matched #f = 11051, uninserted #f = 2128
[2025-02-17 08:58:15.227] [float-tetwild] [info] insert_one_triangle * n done, #v = 14322, #t = 81107
[2025-02-17 08:58:15.227] [float-tetwild] [info] uninserted #f = 0/2128
[2025-02-17 08:58:15.227] [float-tetwild] [info] total timing: 0s
[2025-02-17 08:58:15.248] [float-tetwild] [info] pair_track_surface_fs done
#boundary_e1 = 1
#boundary_e2 = 2
[2025-02-17 08:58:15.252] [float-tetwild] [info] find_boundary_edges done
[2025-02-17 08:58:15.254] [float-tetwild] [info] time1 = 0.002028
[2025-02-17 08:58:15.254] [float-tetwild] [info] uninsert boundary #e = 0/3
[2025-02-17 08:58:15.254] [float-tetwild] [info] time2 = 0
[2025-02-17 08:58:15.254] [float-tetwild] [info] time3 = 5.7e-05
[2025-02-17 08:58:15.254] [float-tetwild] [info] time4 = 0
[2025-02-17 08:58:15.254] [float-tetwild] [info] time5 = 0
[2025-02-17 08:58:15.254] [float-tetwild] [info] time6 = 0
[2025-02-17 08:58:15.254] [float-tetwild] [info] uninserted #f = 0/2128
known_surface_fs.size = 0
known_not_surface_fs.size = 0
[2025-02-17 08:58:15.558] [float-tetwild] [info] mark_surface_fs done
[2025-02-17 08:58:15.559] [float-tetwild] [info] #b_edge1 = 3, #b_edges2 = 0
[2025-02-17 08:58:15.561] [float-tetwild] [info] cutting 0.743994s
[2025-02-17 08:58:15.561] [float-tetwild] [info]
initializing...
edge collapsing...
fixed 0 tangled element
success(env) = 4542
success = 6232(115359)
success(env) = 140
success = 236(78934)
success(env) = 9
success = 65(20400)
success(env) = 0
success = 28(4374)
success(env) = 0
success = 10(1404)
success(env) = 0
success = 3(639)
success(env) = 0
success = 0(115)
edge collapsing done!
time = 0.953942s
#v = 7674
#t = 45793
max_energy = 308182
avg_energy = 18.0234
//////////////// pass 0 ////////////////
edge splitting...
fixed 0 tangled element
success = 740(740)
edge splitting done!
time = 0.024245s
#v = 8414
#t = 49171
max_energy = 308182
avg_energy = 17.0891
edge collapsing...
fixed 0 tangled element
success(env) = 0
success = 534(98185)
success(env) = 0
success = 65(10029)
success(env) = 0
success = 21(3575)
success(env) = 0
success = 8(1027)
success(env) = 0
success = 3(462)
success(env) = 0
success = 1(169)
success(env) = 0
success = 0(9)
edge collapsing done!
time = 0.364032s
#v = 7782
#t = 46321
max_energy = 308182
avg_energy = 17.8367
edge swapping...
fixed 0 tangled element
success3 = 5454
success4 = 5357
success5 = 565
success = 11376(49679)
edge swapping done!
time = 0.250046s
#v = 7782
#t = 41432
max_energy = 429.93
avg_energy = 9.49887
vertex smoothing...
success = 3676(7107)
vertex smoothing done!
time = 0.086589s
#v = 7782
#t = 41432
max_energy = 213.988
avg_energy = 8.32356
//////////////// pass 1 ////////////////
edge splitting...
fixed 0 tangled element
success = 394(394)
edge splitting done!
time = 0.021097s
#v = 8176
#t = 43100
max_energy = 213.988
avg_energy = 8.17282
edge collapsing...
fixed 0 tangled element
success(env) = 433
success = 832(86085)
success(env) = 27
success = 73(40553)
success(env) = 2
success = 14(6655)
success(env) = 0
success = 4(1059)
success(env) = 0
success = 0(152)
edge collapsing done!
time = 0.462798s
#v = 7253
#t = 38778
max_energy = 159.855
avg_energy = 8.16994
edge swapping...
fixed 0 tangled element
success3 = 1012
success4 = 2473
success5 = 303
success = 3788(27374)
edge swapping done!
time = 0.138985s
#v = 7253
#t = 38069
max_energy = 72.6187
avg_energy = 7.84231
vertex smoothing...
success = 2920(6611)
vertex smoothing done!
time = 0.076763s
#v = 7253
#t = 38069
max_energy = 71.3286
avg_energy = 7.59206
//////////////// postprocessing ////////////////
edge collapsing...
fixed 0 tangled element
success(env) = 137
success = 191(79224)
success(env) = 6
success = 15(18112)
success(env) = 3
success = 4(1470)
success(env) = 0
success = 0(360)
edge collapsing done!
time = 0.336444s
#v = 7043
#t = 37065
max_energy = 71.3286
avg_energy = 7.58398
edge collapsing...
fixed 0 tangled element
success(env) = 551
success = 1430(74743)
success(env) = 40
success = 51(49931)
success(env) = 3
success = 3(7971)
success(env) = 0
success = 0(633)
edge collapsing done!
time = 0.731982s
#v = 5559
#t = 30014
max_energy = 99.7259
avg_energy = 10.4634
edge swapping...
fixed 0 tangled element
success3 = 1062
success4 = 2528
success5 = 652
success = 4242(23540)
edge swapping done!
time = 0.140914s
#v = 5559
#t = 29604
max_energy = 91.9226
avg_energy = 9.1182
edge collapsing...
fixed 0 tangled element
success(env) = 95
success = 135(63699)
success(env) = 14
success = 14(17912)
success(env) = 2
success = 2(2128)
success(env) = 0
success = 0(265)
edge collapsing done!
time = 0.486191s
#v = 5408
#t = 28794
max_energy = 91.9226
avg_energy = 9.73409
edge swapping...
fixed 0 tangled element
success3 = 151
success4 = 489
success5 = 162
success = 802(14699)
edge swapping done!
time = 0.069319s
#v = 5408
#t = 28805
max_energy = 91.9226
avg_energy = 9.42783
edge collapsing...
fixed 0 tangled element
success(env) = 11
success = 12(61622)
success(env) = 0
success = 0(2024)
edge collapsing done!
time = 0.328129s
#v = 5396
#t = 28740
max_energy = 91.9226
avg_energy = 9.43665
edge swapping...
fixed 0 tangled element
success3 = 11
success4 = 34
success5 = 4
success = 49(12847)
edge swapping done!
time = 0.052613s
#v = 5396
#t = 28733
max_energy = 91.9226
avg_energy = 9.42744
edge collapsing...
fixed 0 tangled element
success(env) = 0
success = 0(61492)
edge collapsing done!
time = 0.310932s
#v = 5396
#t = 28733
max_energy = 91.9226
avg_energy = 9.42744
edge swapping...
fixed 0 tangled element
success3 = 0
success4 = 0
success5 = 0
success = 0(12719)
edge swapping done!
time = 0.052591s
#v = 5396
#t = 28733
max_energy = 91.9226
avg_energy = 9.42744
[2025-02-17 08:58:20.461] [float-tetwild] [info] mesh optimization 4.89995s
[2025-02-17 08:58:20.461] [float-tetwild] [info]
[2025-02-17 08:58:20.477] [float-tetwild] [info] correct_tracked_surface_orientation done
We want to use this mesh within FEniCSx, so we extract the cells, nodes and volume markers from the mesh instance in wildmeshing.
point_array, cell_array, marker = tetra.get_tet_mesh()
[2025-02-17 08:58:20.582] [float-tetwild] [info] after winding number
[2025-02-17 08:58:20.582] [float-tetwild] [info] #v = 5396
[2025-02-17 08:58:20.582] [float-tetwild] [info] #t = 15117
[2025-02-17 08:58:20.582] [float-tetwild] [info] winding number 1.73978e+09s
[2025-02-17 08:58:20.582] [float-tetwild] [info]
We can pass these into the DOLFINx mesh constructor called create_mesh.
mesh = dolfinx.mesh.create_mesh(
MPI.COMM_WORLD,
cell_array.astype(np.int64),
point_array,
ufl.Mesh(basix.ufl.element("Lagrange", "tetrahedron", 1, shape=(3,))),
)
We add the cell markers as CellTags with the following lines of code
local_entities, local_values = dolfinx.io.gmshio.distribute_entity_data(
mesh,
mesh.topology.dim,
cell_array.astype(np.int64),
marker.flatten().astype(np.int32),
)
adj = dolfinx.graph.adjacencylist(local_entities)
ct = dolfinx.mesh.meshtags_from_entities(
mesh,
mesh.topology.dim,
adj,
local_values.astype(np.int32, copy=False),
)
We visualize the resulting mesh and cell tags with pyvista
pv_grid = pyvista.UnstructuredGrid(*dolfinx.plot.vtk_mesh(mesh))
pv_grid.cell_data["marker"] = ct.values
pyvista.start_xvfb(1.0)
plotter = pyvista.Plotter()
plotter.add_mesh(pv_grid)
plotter.show()
We also create som slices to see the various markers within the mesh
slices = pv_grid.slice_along_axis(n=7, axis="y")
slices.plot()