Generating meshes with Pytetwild

In this demo, we will show how to generate a mesh from STL/ply files based on MRI images. We use the files from the paper In-silico molecular enrichment and clearance of the human intracranial space

We start by importing all the packages required for this task.

import json
import os
from pathlib import Path

from mpi4py import MPI

import numpy as np
import pytetwild
import pyvista

import basix.ufl
import dolfinx
import ufl

Next, we use the already downloaded stl files from this chapter. We will use pytetwild to combine these. We create unique volume markers for the lateral ventricles (LV.ply) and the third and fourth ventricle (V34.ply). Next we embed the ventricular system within the brain parenchyma (parenchyma_incl_ventr.ply), before we add the dura interface (skull.ply), which is the interfacebetween the CSF and the skull.

We extract these files from the "WILDFENICS_DATA_PATH". If you have not set this variable, please change it to the appropriate location on your system.

folder = Path(os.environ["WILDFENICS_DATA_PATH"])
assert folder.exists(), "Could not find surface files"
tree = {
    "operation": "union",
    "left": str((folder / "skull.ply").absolute().as_posix()),
    "right": {
        "operation": "union",
        "left": str((folder / "parenchyma_incl_ventr.ply").absolute().as_posix()),
        "right": {
            "operation": "union",
            "left": str((folder / "LV.ply").absolute().as_posix()),
            "right": str((folder / "V34.ply").absolute().as_posix()),
        },
    },
}

Next, we can set up a 3D meshing instance and load in the instruction-set. We use a named temporary file to pass the options above to pytetwild.

import tempfile

with tempfile.NamedTemporaryFile(mode="w+t", delete=True) as f:
    f.write(json.dumps(tree))
    f.flush()
    tetra = pytetwild.tetrahedralize_csg(
        f.name, stop_energy=10, edge_length_r=0.03, epsilon=0.00225
    )

Pre-checking CSG input meshes for NaN/Inf values...
TBB threads 4
All input meshes passed check.
bbox_diag_length = 0.266221
ideal_edge_length = 0.00798662
stage = 2
eps_input = 0.000598997
eps = 0.000331599
eps_simplification = 0.000265279
eps_coplanar = 2.66221e-07
dd = 0.000399331
dd_simplification = 0.000319465
collapsing 22.1708
swapping 0.293877
#boundary_e1 = 309
#boundary_e2 = 120
FAIL subdivide_tets
#boundary_e1 = 312
#boundary_e2 = 120
known_surface_fs.size = 0
known_not_surface_fs.size = 0
initializing...
edge collapsing...
fixed 0 tangled element
success(env) = 85419
success = 100347(816984)
success(env) = 3247
success = 3461(533465)
success(env) = 284
success = 299(224884)
success(env) = 38
success = 39(28642)
success(env) = 4
success = 4(4807)
success(env) = 1
success = 1(821)
success(env) = 0
success = 0(119)
edge collapsing done!
time = 15.4576s
#v = 44022
#t = 287926
max_energy = 7950.47
avg_energy = 7.75693
#boundary_e1 = 309
#boundary_e2 = 120
known_surface_fs.size = 0
known_not_surface_fs.size = 0
initializing...
//////////////// pass 0 ////////////////
edge spTetrahedralization process from CSG completed successfully.
litting...
fixed 0 tangled element
success = 43033(43033)
edge splitting done!
time = 0.349607s
#v = 87081
#t = 503284
max_energy = 7950.47
avg_energy = 7.68547
edge collapsing...
fixed 0 tangled element
success(env) = 243
success = 29401(628387)
success(env) = 85
success = 1362(165525)
success(env) = 17
success = 189(47746)
success(env) = 2
success = 37(8590)
success(env) = 1
success = 9(1491)
success(env) = 0
success = 2(379)
success(env) = 0
success = 0(109)
edge collapsing done!
time = 5.81219s
#v = 56081
#t = 348660
max_energy = 2689.49
avg_energy = 6.61967
edge swapping...
fixed 0 tangled element
success3 = 33750
success4 = 39777
success5 = 3945
success = 77472(345082)
edge swapping done!
time = 2.7393s
#v = 56081
#t = 318855
max_energy = 1840.36
avg_energy = 5.55799
vertex smoothing...
success = 30042(52182)
vertex smoothing done!
time = 2.03316s
#v = 56081
#t = 318855
max_energy = 514.72
avg_energy = 5.19963
//////////////// pass 1 ////////////////
edge splitting...
fixed 0 tangled element
success = 5071(5071)
edge splitting done!
time = 0.197751s
#v = 61152
#t = 343678
max_energy = 514.72
avg_energy = 5.14713
edge collapsing...
fixed 0 tangled element
success(env) = 2708
success = 6999(524194)
success(env) = 169
success = 378(267513)
success(env) = 16
success = 31(26686)
success(env) = 0
success = 4(2599)
success(env) = 1
success = 4(114)
success(env) = 0
success = 2(153)
success(env) = 0
success = 0(11)
edge collapsing done!
time = 5.55023s
#v = 53734
#t = 305324
max_energy = 71.6353
avg_energy = 5.1033
edge swapping...
fixed 0 tangled element
success3 = 3865
success4 = 13913
success5 = 1189
success = 18967(199343)
edge swapping done!
time = 1.50075s
#v = 53734
#t = 302648
max_energy = 71.6353
avg_energy = 4.99874
vertex smoothing...
success = 25288(49838)
vertex smoothing done!
time = 1.81588s
#v = 53734
#t = 302648
max_energy = 71.6353
avg_energy = 4.88919
//////////////// pass 2 ////////////////
edge splitting...
fixed 0 tangled element
success = 2751(2751)
edge splitting done!
time = 0.17712s
#v = 56485
#t = 315713
max_energy = 71.6353
avg_energy = 4.86681
edge collapsing...
fixed 0 tangled element
success(env) = 875
success = 3393(498334)
success(env) = 53
success = 136(128751)
success(env) = 2
success = 13(9249)
success(env) = 1
success = 2(571)
success(env) = 1
success = 1(170)
success(env) = 0
success = 0(244)
edge collapsing done!
time = 4.07998s
#v = 52940
#t = 298098
max_energy = 71.6353
avg_energy = 4.87446
edge swapping...
fixed 0 tangled element
success3 = 1617
success4 = 7204
success5 = 639
success = 9460(174507)
edge swapping done!
time = 1.1762s
#v = 52940
#t = 297120
max_energy = 71.6353
avg_energy = 4.83822
vertex smoothing...
success = 21984(49044)
vertex smoothing done!
time = 1.65378s
#v = 52940
#t = 297120
max_energy = 71.6353
avg_energy = 4.78795
updating sclaing field ...
filter_energy = 8
is_hit_min_edge_length = 0
enlarge envelope, eps = 0.000368443
//////////////// pass 3 ////////////////
edge splitting...
fixed 0 tangled element
success = 70648(70648)
edge splitting done!
time = 0.518432s
#v = 123588
#t = 674430
max_energy = 71.6353
avg_energy = 4.89067
edge collapsing...
fixed 0 tangled element
success(env) = 2460
success = 47795(555750)
success(env) = 318
success = 3604(309426)
success(env) = 26
success = 417(94643)
success(env) = 5
success = 84(13497)
success(env) = 1
success = 24(2907)
success(env) = 2
success = 5(770)
success(env) = 0
success = 0(218)
edge collapsing done!
time = 7.0936s
#v = 71659
#t = 409647
max_energy = 18.0889
avg_energy = 4.37828
edge swapping...
fixed 0 tangled element
success3 = 5361
success4 = 23730
success5 = 2124
success = 31215(305017)
edge swapping done!
time = 2.92811s
#v = 71659
#t = 406410
max_energy = 18.0889
avg_energy = 4.21998
vertex smoothing...
success = 46040(67748)
vertex smoothing done!
time = 2.15115s
#v = 71659
#t = 406410
max_energy = 17.4131
avg_energy = 4.02804
//////////////// pass 4 ////////////////
edge splitting...
fixed 0 tangled element
success = 15201(15201)
edge splitting done!
time = 0.28983s
#v = 86860
#t = 480330
max_energy = 19.2039
avg_energy = 4.06878
edge collapsing...
fixed 0 tangled element
success(env) = 1415
success = 14825(435219)
success(env) = 94
success = 688(192714)
success(env) = 8
success = 68(19700)
success(env) = 0
success = 4(1858)
success(env) = 0
success = 2(62)
success(env) = 0
success = 0(9)
edge collapsing done!
time = 4.24873s
#v = 71273
#t = 404532
max_energy = 17.4131
avg_energy = 4.02155
edge swapping...
fixed 0 tangled element
success3 = 1530
success4 = 11818
success5 = 777
success = 14125(269261)
edge swapping done!
time = 1.92521s
#v = 71273
#t = 403779
max_energy = 17.4131
avg_energy = 3.98579
vertex smoothing...
success = 43328(67355)
vertex smoothing done!
time = 2.03905s
#v = 71273
#t = 403779
max_energy = 17.4131
avg_energy = 3.90952
updating sclaing field ...
filter_energy = 8
is_hit_min_edge_length = 0
//////////////// pass 5 ////////////////
edge splitting...
fixed 0 tangled element
success = 27156(27156)
edge splitting done!
time = 0.351325s
#v = 98429
#t = 544773
max_energy = 17.6169
avg_energy = 4.06145
edge collapsing...
fixed 0 tangled element
success(env) = 2499
success = 18783(637947)
success(env) = 196
success = 1433(203831)
success(env) = 15
success = 175(27951)
success(env) = 3
success = 39(4051)
success(env) = 1
success = 7(922)
success(env) = 0
success = 1(162)
success(env) = 0
success = 0(10)
edge collapsing done!
time = 5.30941s
#v = 77991
#t = 441978
max_energy = 10.1712
avg_energy = 3.91513
edge swapping...
fixed 0 tangled element
success3 = 1258
success4 = 10288
success5 = 684
success = 12230(294192)
edge swapping done!
time = 2.31537s
#v = 77991
#t = 441404
max_energy = 9.5359
avg_energy = 3.88536
vertex smoothing...
success = 48703(74071)
vertex smoothing done!
time = 2.09602s
#v = 77991
#t = 441404
max_energy = 9.45595
avg_energy = 3.81207
//////////////// postprocessing ////////////////
edge collapsing...
fixed 0 tangled element
success(env) = 421
success = 1727(882652)
success(env) = 25
success = 143(154157)
success(env) = 4
success = 30(15395)
success(env) = 1
success = 6(3372)
success(env) = 0
success = 1(642)
success(env) = 0
success = 0(219)
edge collapsing done!
time = 5.38895s
#v = 76084
#t = 430997
max_energy = 9.45595
avg_energy = 3.81348
edge collapsing...
fixed 0 tangled element
success(env) = 6547
success = 28997(572474)
success(env) = 445
success = 910(424804)
success(env) = 38
success = 89(95891)
success(env) = 4
success = 9(10348)
success(env) = 0
success = 0(1215)
edge collapsing done!
time = 7.75637s
#v = 46079
#t = 259317
max_energy = 9.99985
avg_energy = 4.81062
edge swapping...
fixed 0 tangled element
success3 = 5314
success4 = 19600
success5 = 2129
success = 27043(189106)
edge swapping done!
time = 1.27421s
#v = 46079
#t = 256132
max_energy = 9.99919
avg_energy = 4.56349
edge collapsing...
fixed 0 tangled element
success(env) = 839
success = 2381(412552)
success(env) = 70
success = 129(195529)
success(env) = 7
success = 11(16444)
success(env) = 1
success = 1(1287)
success(env) = 0
success = 0(119)
edge collapsing done!
time = 4.01351s
#v = 43557
#t = 241038
max_energy = 9.99989
avg_energy = 4.7272
edge swapping...
fixed 0 tangled element
success3 = 1164
success4 = 4794
success5 = 710
success = 6668(135889)
edge swapping done!
time = 0.84189s
#v = 43557
#t = 240584
max_energy = 9.99823
avg_energy = 4.66804
edge collapsing...
fixed 0 tangled element
success(env) = 94
success = 353(389817)
success(env) = 11
success = 22(42235)
success(env) = 1
success = 3(3271)
success(env) = 2
success = 2(474)
success(env) = 0
success = 0(260)
edge collapsing done!
time = 2.74564s
#v = 43177
#t = 238242
max_energy = 9.99888
avg_energy = 4.69602
edge swapping...
fixed 0 tangled element
success3 = 183
success4 = 879
success5 = 139
success = 1201(122479)
edge swapping done!
time = 0.715619s
#v = 43177
#t = 238198
max_energy = 9.99888
avg_energy = 4.68595
edge collapsing...
fixed 0 tangled element
success(env) = 5
success = 34(386238)
success(env) = 1
success = 3(4978)
success(env) = 0
success = 0(325)
edge collapsing done!
time = 2.51101s
#v = 43140
#t = 237974
max_energy = 9.99888
avg_energy = 4.68925
edge swapping...
fixed 0 tangled element
success3 = 28
success4 = 138
success5 = 20
success = 186(120147)
edge swapping done!
time = 0.706955s
#v = 43140
#t = 237966
max_energy = 9.99888
avg_energy = 4.68762
optimization finished
correct surface orientation finished
boolean operation finished 

We want to use this mesh within FEniCSx, so we extract the cells, nodes and volume markers from tetra, a pyvista.UnstructuredGrid.

cell_array = tetra.cells.reshape(
    -1, 5
)[
    :, 1:
]  # Pyvista store (N0, node1, ..., nodeN0, N1, node1, ..., nodeN1, ...) so we need to skip the first column
point_array = tetra.points.reshape(-1, 3)
marker = tetra.cell_data["marker"]

We can pass these into the DOLFINx mesh constructor called create_mesh.

mesh = dolfinx.mesh.create_mesh(
    MPI.COMM_WORLD,
    cells=cell_array.astype(np.int64),
    x=point_array,
    e=ufl.Mesh(basix.ufl.element("Lagrange", "tetrahedron", 1, shape=(3,))),
)

We add the cell markers as MeshTags with the following lines of code

local_entities, local_values = dolfinx.io.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 store the mesh and cell markers to file for usage in later tutorials

with dolfinx.io.XDMFFile(mesh.comm, folder / "brain.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_meshtags(ct, mesh.geometry)

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

plotter = pyvista.Plotter()
plotter.add_mesh(pv_grid)
plotter.show()

2026-08-14 09:04:33.839 ( 143.221s) [    7F1548066140]vtkXOpenGLRenderWindow.:1460  WARN| bad X server connection. DISPLAY=

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()