Generating meshes with Wildmeshing

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.

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

tetra = wildmeshing.Tetrahedralizer(
    stop_quality=10, max_its=30, edge_length_r=0.03, epsilon=0.00225
)
tetra.load_csg_tree(json.dumps(tree))

TBB threads 4
o-[I/O         ] Loading file /root/data/mesh/standard/surfaces/skull.ply...
                 (FP64) nb_v:138154 nb_e:0 nb_f:276304 nb_b:0 tri:1 dim:3
                 Attributes on vertices: point[3]
                 Loading file /root/data/mesh/standard/surfaces/parenchyma_incl_ventr.ply...
                 (FP64) nb_v:317019 nb_e:0 nb_f:635594 nb_b:0 tri:1 dim:3
                 Attributes on vertices: point[3]
                 Loading file /root/data/mesh/standard/surfaces/LV.ply...
                 (FP64) nb_v:17572 nb_e:0 nb_f:35144 nb_b:0 tri:1 dim:3
                 Attributes on vertices: point[3]
                 Loading file /root/data/mesh/standard/surfaces/V34.ply...
                 (FP64) nb_v:4608 nb_e:0 nb_f:9212 nb_b:0 tri:1 dim:3
                 Attributes on vertices: point[3]
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

Then, we can generate the mesh

# Create mesh
tetra.tetrahedralize()

[2025-10-29 10:04:25.834] [float-tetwild] [info] remove duplicates: 
[2025-10-29 10:04:25.834] [float-tetwild] [info] #v: 477353 -> 477353
[2025-10-29 10:04:25.834] [float-tetwild] [info] #f: 956254 -> 956254
collapsing 21.0959
swapping 0.261264
[2025-10-29 10:04:47.431] [float-tetwild] [info] remove duplicates: 
[2025-10-29 10:04:47.431] [float-tetwild] [info] #v: 49442 -> 49442
[2025-10-29 10:04:47.431] [float-tetwild] [info] #f: 100369 -> 100369
[2025-10-29 10:04:47.433] [float-tetwild] [info] #v = 49442
[2025-10-29 10:04:47.433] [float-tetwild] [info] #f = 100256
#boundary_e1 = 309
#boundary_e2 = 120
[2025-10-29 10:04:47.797] [float-tetwild] [info] preprocessing 22.1805s
[2025-10-29 10:04:47.797] [float-tetwild] [info] 
[2025-10-29 10:04:48.029] [float-tetwild] [info] #v = 50230
[2025-10-29 10:04:48.030] [float-tetwild] [info] #t = 333917
[2025-10-29 10:04:48.030] [float-tetwild] [info] tetrahedralizing 0.232586s
[2025-10-29 10:04:48.030] [float-tetwild] [info] 
[2025-10-29 10:04:48.031] [float-tetwild] [info] triangle insertion start, #f = 100256, #v = 50230, #t = 333917
[2025-10-29 10:04:48.166] [float-tetwild] [info] matched #f = 74494, uninserted #f = 25762
[2025-10-29 10:04:57.068] [float-tetwild] [info] insert_one_triangle * n done, #v = 153581, #t = 867003
[2025-10-29 10:04:57.068] [float-tetwild] [info] uninserted #f = 0/25762
[2025-10-29 10:04:57.068] [float-tetwild] [info] total timing: 0s
[2025-10-29 10:04:57.510] [float-tetwild] [info] pair_track_surface_fs done
#boundary_e1 = 309
#boundary_e2 = 120
[2025-10-29 10:04:57.560] [float-tetwild] [info] find_boundary_edges done
[2025-10-29 10:04:57.608] [float-tetwild] [info] time1 = 0.047971
[2025-10-29 10:04:57.620] [float-tetwild] [info] uninsert boundary #e = 0/429
[2025-10-29 10:04:57.620] [float-tetwild] [info] time2 = 4.7e-05
[2025-10-29 10:04:57.620] [float-tetwild] [info] time3 = 0.011125
[2025-10-29 10:04:57.620] [float-tetwild] [info] time4 = 0.000278
[2025-10-29 10:04:57.620] [float-tetwild] [info] time5 = 7.2e-05
[2025-10-29 10:04:57.620] [float-tetwild] [info] time6 = 1e-06
[2025-10-29 10:04:57.626] [float-tetwild] [info] uninserted #f = 0/25762
known_surface_fs.size = 0
known_not_surface_fs.size = 0
[2025-10-29 10:05:00.307] [float-tetwild] [info] mark_surface_fs done
[2025-10-29 10:05:00.329] [float-tetwild] [info] #b_edge1 = 429, #b_edges2 = 0
[2025-10-29 10:05:00.361] [float-tetwild] [info] cutting 12.3301s
[2025-10-29 10:05:00.361] [float-tetwild] [info] 
initializing...
edge collapsing...
fixed 0 tangled element
success(env) = 85287
success = 100284(815086)
success(env) = 3254
success = 3481(532493)
success(env) = 258
success = 267(232191)
success(env) = 41
success = 42(27456)
success(env) = 7
success = 7(5093)
success(env) = 3
success = 4(1205)
success(env) = 0
success = 0(711)
edge collapsing done!
time = 14.0169s
#v = 43980
#t = 287645
max_energy = 7703.34
avg_energy = 7.7571
//////////////// pass 0 ////////////////
edge splitting...
fixed 0 tangled element
success = 43060(43060)
edge splitting done!
time = 0.324342s
#v = 87040
#t = 503038
max_energy = 7703.34
avg_energy = 7.71727
edge collapsing...
fixed 0 tangled element
success(env) = 180
success = 29386(621377)
success(env) = 88
success = 1404(157444)
success(env) = 25
success = 207(48741)
success(env) = 8
success = 44(9129)
success(env) = 1
success = 18(2131)
success(env) = 0
success = 3(823)
success(env) = 0
success = 0(122)
edge collapsing done!
time = 5.27286s
#v = 55978
#t = 348053
max_energy = 3132.41
avg_energy = 6.63385
edge swapping...
fixed 0 tangled element
success3 = 33769
success4 = 39409
success5 = 3871
success = 77049(343192)
edge swapping done!
time = 2.54954s
#v = 55978
#t = 318155
max_energy = 1749.69
avg_energy = 5.55675
vertex smoothing...
success = 29863(52083)
vertex smoothing done!
time = 1.99577s
#v = 55978
#t = 318155
max_energy = 496.105
avg_energy = 5.21014
//////////////// pass 1 ////////////////
edge splitting...
fixed 0 tangled element
success = 5200(5200)
edge splitting done!
time = 0.197336s
#v = 61178
#t = 343624
max_energy = 496.105
avg_energy = 5.15865
edge collapsing...
fixed 0 tangled element
success(env) = 2654
success = 6886(518932)
success(env) = 149
success = 332(260480)
success(env) = 11
success = 28(22487)
success(env) = 0
success = 2(1737)
success(env) = 0
success = 1(34)
success(env) = 0
success = 0(20)
edge collapsing done!
time = 4.6982s
#v = 53929
#t = 306335
max_energy = 245.189
avg_energy = 5.11177
edge swapping...
fixed 0 tangled element
success3 = 3915
success4 = 14106
success5 = 1193
success = 19214(200571)
edge swapping done!
time = 1.36475s
#v = 53929
#t = 303613
max_energy = 245.189
avg_energy = 5.00667
vertex smoothing...
success = 25321(50035)
vertex smoothing done!
time = 1.75813s
#v = 53929
#t = 303613
max_energy = 244.251
avg_energy = 4.89625
//////////////// pass 2 ////////////////
edge splitting...
fixed 0 tangled element
success = 2678(2678)
edge splitting done!
time = 0.174516s
#v = 56607
#t = 316327
max_energy = 244.251
avg_energy = 4.87565
edge collapsing...
fixed 0 tangled element
success(env) = 847
success = 3258(494608)
success(env) = 43
success = 99(122070)
success(env) = 3
success = 11(6325)
success(env) = 2
success = 2(545)
success(env) = 0
success = 0(305)
edge collapsing done!
time = 3.42877s
#v = 53237
#t = 299679
max_energy = 62.1568
avg_energy = 4.87551
edge swapping...
fixed 0 tangled element
success3 = 1547
success4 = 6998
success5 = 551
success = 9096(174015)
edge swapping done!
time = 1.08824s
#v = 53237
#t = 298683
max_energy = 62.1568
avg_energy = 4.83906
vertex smoothing...
success = 21791(49342)
vertex smoothing done!
time = 1.61941s
#v = 53237
#t = 298683
max_energy = 61.6505
avg_energy = 4.79242
//////////////// pass 3 ////////////////
edge splitting...
fixed 0 tangled element
success = 1903(1903)
edge splitting done!
time = 0.168367s
#v = 55140
#t = 307580
max_energy = 61.6505
avg_energy = 4.77872
edge collapsing...
fixed 0 tangled element
success(env) = 390
success = 2208(486481)
success(env) = 22
success = 59(73837)
success(env) = 4
success = 6(4192)
success(env) = 0
success = 1(474)
success(env) = 0
success = 0(31)
edge collapsing done!
time = 3.0508s
#v = 52866
#t = 296637
max_energy = 61.6505
avg_energy = 4.78848
edge swapping...
fixed 0 tangled element
success3 = 819
success4 = 3946
success5 = 311
success = 5076(163710)
edge swapping done!
time = 0.993939s
#v = 52866
#t = 296129
max_energy = 61.6505
avg_energy = 4.77204
vertex smoothing...
success = 19411(48967)
vertex smoothing done!
time = 1.52327s
#v = 52866
#t = 296129
max_energy = 61.6505
avg_energy = 4.74717
updating sclaing field ...
filter_energy = 8
is_hit_min_edge_length = 0
enlarge envelope, eps = 0.000368443
//////////////// pass 4 ////////////////
edge splitting...
fixed 0 tangled element
success = 66244(66244)
edge splitting done!
time = 0.470816s
#v = 119110
#t = 649793
max_energy = 55.0953
avg_energy = 4.86133
edge collapsing...
fixed 0 tangled element
success(env) = 2250
success = 45076(537159)
success(env) = 271
success = 3338(295709)
success(env) = 41
success = 437(86721)
success(env) = 13
success = 85(14643)
success(env) = 2
success = 16(2978)
success(env) = 0
success = 4(696)
success(env) = 0
success = 0(58)
edge collapsing done!
time = 6.11582s
#v = 70154
#t = 400119
max_energy = 55.0953
avg_energy = 4.3707
edge swapping...
fixed 0 tangled element
success3 = 4905
success4 = 21462
success5 = 1986
success = 28353(291967)
edge swapping done!
time = 2.63282s
#v = 70154
#t = 397200
max_energy = 55.0953
avg_energy = 4.22386
vertex smoothing...
success = 43846(66239)
vertex smoothing done!
time = 2.04963s
#v = 70154
#t = 397200
max_energy = 21.6008
avg_energy = 4.0436
//////////////// pass 5 ////////////////
edge splitting...
fixed 0 tangled element
success = 14713(14713)
edge splitting done!
time = 0.277941s
#v = 84867
#t = 468664
max_energy = 24.0827
avg_energy = 4.08235
edge collapsing...
fixed 0 tangled element
success(env) = 1341
success = 14475(425723)
success(env) = 79
success = 607(183937)
success(env) = 9
success = 47(17506)
success(env) = 0
success = 3(1518)
success(env) = 0
success = 0(64)
edge collapsing done!
time = 3.64827s
#v = 69735
#t = 395112
max_energy = 21.6008
avg_energy = 4.03901
edge swapping...
fixed 0 tangled element
success3 = 1458
success4 = 10888
success5 = 773
success = 13119(259627)
edge swapping done!
time = 1.7323s
#v = 69735
#t = 394427
max_energy = 21.6008
avg_energy = 4.00503
vertex smoothing...
success = 41483(65820)
vertex smoothing done!
time = 1.92295s
#v = 69735
#t = 394427
max_energy = 21.1336
avg_energy = 3.93325
updating sclaing field ...
filter_energy = 8
is_hit_min_edge_length = 0
//////////////// pass 6 ////////////////
edge splitting...
fixed 0 tangled element
success = 34445(34445)
edge splitting done!
time = 0.36129s
#v = 104180
#t = 574157
max_energy = 26.5526
avg_energy = 4.10965
edge collapsing...
fixed 0 tangled element
success(env) = 2834
success = 22850(619631)
success(env) = 231
success = 1690(219547)
success(env) = 25
success = 209(33474)
success(env) = 1
success = 45(4984)
success(env) = 1
success = 8(1384)
success(env) = 1
success = 3(236)
success(env) = 1
success = 1(92)
success(env) = 0
success = 0(22)
edge collapsing done!
time = 5.07088s
#v = 79374
#t = 449324
max_energy = 16.3974
avg_energy = 3.94028
edge swapping...
fixed 0 tangled element
success3 = 1440
success4 = 11396
success5 = 735
success = 13571(302506)
edge swapping done!
time = 2.46078s
#v = 79374
#t = 448619
max_energy = 16.3974
avg_energy = 3.90466
vertex smoothing...
success = 49793(75464)
vertex smoothing done!
time = 2.07369s
#v = 79374
#t = 448619
max_energy = 15.3876
avg_energy = 3.81909
//////////////// pass 7 ////////////////
edge splitting...
fixed 0 tangled element
success = 11665(11665)
edge splitting done!
time = 0.290284s
#v = 91039
#t = 505558
max_energy = 15.3876
avg_energy = 3.87836
edge collapsing...
fixed 0 tangled element
success(env) = 1588
success = 11636(563325)
success(env) = 71
success = 477(141316)
success(env) = 8
success = 58(11516)
success(env) = 1
success = 6(1425)
success(env) = 0
success = 0(79)
edge collapsing done!
time = 3.74452s
#v = 78862
#t = 446323
max_energy = 15.3876
avg_energy = 3.8239
edge swapping...
fixed 0 tangled element
success3 = 529
success4 = 6227
success5 = 324
success = 7080(288567)
edge swapping done!
time = 1.77254s
#v = 78862
#t = 446118
max_energy = 15.3876
avg_energy = 3.81397
vertex smoothing...
success = 47329(74947)
vertex smoothing done!
time = 1.98324s
#v = 78862
#t = 446118
max_energy = 15.3371
avg_energy = 3.77634
updating sclaing field ...
filter_energy = 8
is_hit_min_edge_length = 0
//////////////// pass 8 ////////////////
edge splitting...
fixed 0 tangled element
success = 2475(2475)
edge splitting done!
time = 0.248976s
#v = 81337
#t = 458742
max_energy = 17.1503
avg_energy = 3.79374
edge collapsing...
fixed 0 tangled element
success(env) = 454
success = 2638(823809)
success(env) = 38
success = 178(92152)
success(env) = 5
success = 22(8162)
success(env) = 0
success = 10(1238)
success(env) = 0
success = 5(428)
success(env) = 1
success = 2(215)
success(env) = 0
success = 1(65)
success(env) = 0
success = 0(47)
edge collapsing done!
time = 4.78439s
#v = 78481
#t = 444137
max_energy = 10.3452
avg_energy = 3.78161
edge swapping...
fixed 0 tangled element
success3 = 219
success4 = 3250
success5 = 111
success = 3580(280925)
edge swapping done!
time = 1.66185s
#v = 78481
#t = 444029
max_energy = 10.3452
avg_energy = 3.77776
vertex smoothing...
success = 45100(74564)
vertex smoothing done!
time = 1.90859s
#v = 78481
#t = 444029
max_energy = 10.1918
avg_energy = 3.7599
//////////////// pass 9 ////////////////
edge splitting...
fixed 0 tangled element
success = 1323(1323)
edge splitting done!
time = 0.24593s
#v = 79804
#t = 450352
max_energy = 10.1918
avg_energy = 3.76744
edge collapsing...
fixed 0 tangled element
success(env) = 252
success = 1665(814274)
success(env) = 12
success = 61(40978)
success(env) = 2
success = 8(2683)
success(env) = 0
success = 0(468)
edge collapsing done!
time = 4.07154s
#v = 78070
#t = 441809
max_energy = 10.1918
avg_energy = 3.76167
edge swapping...
fixed 0 tangled element
success3 = 117
success4 = 1679
success5 = 64
success = 1860(276198)
edge swapping done!
time = 1.5684s
#v = 78070
#t = 441756
max_energy = 10.1918
avg_energy = 3.76024
vertex smoothing...
success = 43147(74162)
vertex smoothing done!
time = 1.83065s
#v = 78070
#t = 441756
max_energy = 9.91435
avg_energy = 3.75201
updating sclaing field ...
filter_energy = 8
is_hit_min_edge_length = 1
//////////////// postprocessing ////////////////
edge collapsing...
fixed 0 tangled element
success(env) = 49
success = 202(890449)
success(env) = 4
success = 8(21775)
success(env) = 1
success = 2(865)
success(env) = 0
success = 0(144)
edge collapsing done!
time = 4.24761s
#v = 77858
#t = 440586
max_energy = 8.99006
avg_energy = 3.75246
edge collapsing...
fixed 0 tangled element
success(env) = 6773
success = 30911(573201)
success(env) = 471
success = 939(418815)
success(env) = 28
success = 56(94980)
success(env) = 4
success = 5(6833)
success(env) = 2
success = 2(638)
success(env) = 0
success = 0(290)
edge collapsing done!
time = 7.78731s
#v = 45945
#t = 258197
max_energy = 9.99994
avg_energy = 4.80482
edge swapping...
fixed 0 tangled element
success3 = 5188
success4 = 19562
success5 = 2151
success = 26901(188416)
edge swapping done!
time = 1.26615s
#v = 45945
#t = 255160
max_energy = 9.99976
avg_energy = 4.55948
edge collapsing...
fixed 0 tangled element
success(env) = 842
success = 2435(405592)
success(env) = 62
success = 116(194980)
success(env) = 3
success = 4(13955)
success(env) = 0
success = 1(444)
success(env) = 0
success = 0(284)
edge collapsing done!
time = 3.90521s
#v = 43389
#t = 239794
max_energy = 9.99976
avg_energy = 4.72874
edge swapping...
fixed 0 tangled element
success3 = 1162
success4 = 4684
success5 = 712
success = 6558(134803)
edge swapping done!
time = 0.840007s
#v = 43389
#t = 239344
max_energy = 9.99976
avg_energy = 4.67002
edge collapsing...
fixed 0 tangled element
success(env) = 100
success = 360(382414)
success(env) = 6
success = 19(44041)
success(env) = 2
success = 3(2999)
success(env) = 0
success = 0(458)
edge collapsing done!
time = 2.70525s
#v = 43007
#t = 236994
max_energy = 9.99976
avg_energy = 4.70163
edge swapping...
fixed 0 tangled element
success3 = 219
success4 = 993
success5 = 163
success = 1375(121981)
edge swapping done!
time = 0.714702s
#v = 43007
#t = 236938
max_energy = 9.99976
avg_energy = 4.68981
edge collapsing...
fixed 0 tangled element
success(env) = 8
success = 49(378614)
success(env) = 0
success = 2(6650)
success(env) = 1
success = 2(311)
success(env) = 0
success = 0(194)
edge collapsing done!
time = 2.45007s
#v = 42954
#t = 236601
max_energy = 9.99976
avg_energy = 4.69456
edge swapping...
fixed 0 tangled element
success3 = 16
success4 = 153
success5 = 24
success = 193(119199)
edge swapping done!
time = 0.682352s
#v = 42954
#t = 236609
max_energy = 9.99976
avg_energy = 4.69305
[2025-10-29 10:07:02.734] [float-tetwild] [info] mesh optimization 122.365s
[2025-10-29 10:07:02.734] [float-tetwild] [info] 
[2025-10-29 10:07:02.969] [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-10-29 10:07:04.271] [float-tetwild] [info] after winding number
[2025-10-29 10:07:04.272] [float-tetwild] [info] #v = 42954
[2025-10-29 10:07:04.274] [float-tetwild] [info] #t = 185279
[2025-10-29 10:07:04.274] [float-tetwild] [info] winding number 1.76173e+09s
[2025-10-29 10:07:04.274] [float-tetwild] [info] 

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

2025-10-29 10:07:05.555 ( 161.179s) [    7F4DB4572140]vtkXOpenGLRenderWindow.:1458  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()