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-05-05 13:09:39.722] [float-tetwild] [info] remove duplicates: 
[2025-05-05 13:09:39.722] [float-tetwild] [info] #v: 477353 -> 477353
[2025-05-05 13:09:39.722] [float-tetwild] [info] #f: 956254 -> 956254
collapsing 20.9883
swapping 0.236966
[2025-05-05 13:10:01.182] [float-tetwild] [info] remove duplicates: 
[2025-05-05 13:10:01.182] [float-tetwild] [info] #v: 49442 -> 49442
[2025-05-05 13:10:01.182] [float-tetwild] [info] #f: 100369 -> 100369
[2025-05-05 13:10:01.183] [float-tetwild] [info] #v = 49442
[2025-05-05 13:10:01.183] [float-tetwild] [info] #f = 100256
#boundary_e1 = 309
#boundary_e2 = 120
[2025-05-05 13:10:01.536] [float-tetwild] [info] preprocessing 22.0308s
[2025-05-05 13:10:01.536] [float-tetwild] [info] 
[2025-05-05 13:10:01.767] [float-tetwild] [info] #v = 50230
[2025-05-05 13:10:01.768] [float-tetwild] [info] #t = 333917
[2025-05-05 13:10:01.768] [float-tetwild] [info] tetrahedralizing 0.231984s
[2025-05-05 13:10:01.768] [float-tetwild] [info] 
[2025-05-05 13:10:01.769] [float-tetwild] [info] triangle insertion start, #f = 100256, #v = 50230, #t = 333917
[2025-05-05 13:10:01.904] [float-tetwild] [info] matched #f = 74494, uninserted #f = 25762
[2025-05-05 13:10:10.825] [float-tetwild] [info] insert_one_triangle * n done, #v = 154013, #t = 869356
[2025-05-05 13:10:10.825] [float-tetwild] [info] uninserted #f = 0/25762
[2025-05-05 13:10:10.825] [float-tetwild] [info] total timing: 0s
[2025-05-05 13:10:11.257] [float-tetwild] [info] pair_track_surface_fs done
#boundary_e1 = 309
#boundary_e2 = 120
[2025-05-05 13:10:11.306] [float-tetwild] [info] find_boundary_edges done
[2025-05-05 13:10:11.354] [float-tetwild] [info] time1 = 0.048015
[2025-05-05 13:10:11.366] [float-tetwild] [info] uninsert boundary #e = 0/429
[2025-05-05 13:10:11.366] [float-tetwild] [info] time2 = 4.8e-05
[2025-05-05 13:10:11.366] [float-tetwild] [info] time3 = 0.011399
[2025-05-05 13:10:11.366] [float-tetwild] [info] time4 = 0.000273
[2025-05-05 13:10:11.366] [float-tetwild] [info] time5 = 6.8e-05
[2025-05-05 13:10:11.366] [float-tetwild] [info] time6 = 0
[2025-05-05 13:10:11.372] [float-tetwild] [info] uninserted #f = 0/25762
known_surface_fs.size = 0
known_not_surface_fs.size = 0
[2025-05-05 13:10:14.039] [float-tetwild] [info] mark_surface_fs done
[2025-05-05 13:10:14.061] [float-tetwild] [info] #b_edge1 = 429, #b_edges2 = 0
[2025-05-05 13:10:14.092] [float-tetwild] [info] cutting 12.3236s
[2025-05-05 13:10:14.092] [float-tetwild] [info] 
initializing...
edge collapsing...
fixed 0 tangled element
success(env) = 85383
success = 100384(816175)
success(env) = 3364
success = 3588(532337)
success(env) = 244
success = 252(226250)
success(env) = 37
success = 39(28344)
success(env) = 7
success = 7(4619)
success(env) = 0
success = 0(684)
edge collapsing done!
time = 13.6417s
#v = 43976
#t = 287615
max_energy = 1299.55
avg_energy = 7.63885
//////////////// pass 0 ////////////////
edge splitting...
fixed 0 tangled element
success = 42909(42909)
edge splitting done!
time = 0.320129s
#v = 86885
#t = 502195
max_energy = 1299.55
avg_energy = 7.55227
edge collapsing...
fixed 0 tangled element
success(env) = 173
success = 29186(621095)
success(env) = 78
success = 1401(156965)
success(env) = 18
success = 217(47637)
success(env) = 3
success = 45(7879)
success(env) = 1
success = 16(2032)
success(env) = 0
success = 3(656)
success(env) = 0
success = 0(23)
edge collapsing done!
time = 4.95269s
#v = 56017
#t = 348312
max_energy = 1299.55
avg_energy = 6.60196
edge swapping...
fixed 0 tangled element
success3 = 33739
success4 = 39555
success5 = 3800
success = 77094(343488)
edge swapping done!
time = 2.37469s
#v = 56017
#t = 318373
max_energy = 892.941
avg_energy = 5.55161
vertex smoothing...
success = 30122(52121)
vertex smoothing done!
time = 1.99584s
#v = 56017
#t = 318373
max_energy = 892.941
avg_energy = 5.21507
//////////////// pass 1 ////////////////
edge splitting...
fixed 0 tangled element
success = 5233(5233)
edge splitting done!
time = 0.201772s
#v = 61250
#t = 343846
max_energy = 892.941
avg_energy = 5.16069
edge collapsing...
fixed 0 tangled element
success(env) = 2633
success = 6952(519441)
success(env) = 159
success = 345(261950)
success(env) = 16
success = 41(24993)
success(env) = 3
success = 7(2582)
success(env) = 0
success = 1(334)
success(env) = 0
success = 0(11)
edge collapsing done!
time = 4.45704s
#v = 53904
#t = 306175
max_energy = 126.469
avg_energy = 5.11403
edge swapping...
fixed 0 tangled element
success3 = 3808
success4 = 13847
success5 = 1236
success = 18891(199439)
edge swapping done!
time = 1.3118s
#v = 53904
#t = 303603
max_energy = 126.469
avg_energy = 5.01214
vertex smoothing...
success = 25322(50010)
vertex smoothing done!
time = 1.75277s
#v = 53904
#t = 303603
max_energy = 121.648
avg_energy = 4.90469
//////////////// pass 2 ////////////////
edge splitting...
fixed 0 tangled element
success = 2734(2734)
edge splitting done!
time = 0.181501s
#v = 56638
#t = 316605
max_energy = 121.648
avg_energy = 4.88247
edge collapsing...
fixed 0 tangled element
success(env) = 882
success = 3426(493986)
success(env) = 52
success = 118(128401)
success(env) = 6
success = 14(8222)
success(env) = 0
success = 0(754)
edge collapsing done!
time = 3.33623s
#v = 53080
#t = 298907
max_energy = 121.648
avg_energy = 4.88794
edge swapping...
fixed 0 tangled element
success3 = 1640
success4 = 7186
success5 = 609
success = 9435(174744)
edge swapping done!
time = 1.06938s
#v = 53080
#t = 297876
max_energy = 121.648
avg_energy = 4.85012
vertex smoothing...
success = 21767(49179)
vertex smoothing done!
time = 1.62504s
#v = 53080
#t = 297876
max_energy = 117.391
avg_energy = 4.80133
//////////////// pass 3 ////////////////
edge splitting...
fixed 0 tangled element
success = 2023(2023)
edge splitting done!
time = 0.176111s
#v = 55103
#t = 307356
max_energy = 117.391
avg_energy = 4.78725
edge collapsing...
fixed 0 tangled element
success(env) = 388
success = 2291(485239)
success(env) = 18
success = 60(75595)
success(env) = 2
success = 9(3087)
success(env) = 0
success = 4(446)
success(env) = 0
success = 1(45)
success(env) = 0
success = 0(6)
edge collapsing done!
time = 3.0514s
#v = 52738
#t = 295974
max_energy = 117.391
avg_energ

y = 4.79637
edge swapping...
fixed 0 tangled element
success3 = 803
success4 = 4049
success5 = 320
success = 5172(163565)
edge swapping done!
time = 0.964905s
#v = 52738
#t = 295491
max_energy = 117.391
avg_energy = 4.7793
vertex smoothing...
success = 19411(48841)
vertex smoothing done!
time = 1.52413s
#v = 52738
#t = 295491
max_energy = 117.391
avg_energy = 4.75371
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 = 67125(67125)
edge splitting done!
time = 0.466925s
#v = 119863
#t = 654057
max_energy = 103.261
avg_energy = 4.8618
edge collapsing...
fixed 0 tangled element
success(env) = 2334
success = 45397(540006)
success(env) = 249
success = 3469(296235)
success(env) = 30
success = 368(87368)
success(env) = 7
success = 93(12141)
success(env) = 1
success = 11(2970)
success(env) = 0
success = 2(319)
success(env) = 0
success = 0(73)
edge collapsing done!
time = 6.06135s
#v = 70523
#t = 402571
max_energy = 103.261
avg_energy = 4.36898
edge swapping...
fixed 0 tangled element
success3 = 4995
success4 = 21363
success5 = 1949
success = 28307(293772)
edge swapping done!
time = 2.61262s
#v = 70523
#t = 399525
max_energy = 103.261
avg_energy = 4.2227
vertex smoothing...
success = 44021(66616)
vertex smoothing done!
time = 2.04834s
#v = 70523
#t = 399525
max_energy = 103.261
avg_energy = 4.04145
//////////////// pass 5 ////////////////
edge splitting...
fixed 0 tangled element
success = 14788(14788)
edge splitting done!
time = 0.288109s
#v = 85311
#t = 471341
max_energy = 103.261
avg_energy = 4.0809
edge collapsing...
fixed 0 tangled element
success(env) = 1422
success = 14578(426409)
success(env) = 71
success = 579(186020)
success(env) = 5
success = 43(16009)
success(env) = 1
success = 5(1374)
success(env) = 0
success = 0(228)
edge collapsing done!
time = 3.57971s
#v = 70106
#t = 397469
max_energy = 103.261
avg_energy = 4.03679
edge swapping...
fixed 0 tangled element
success3 = 1488
success4 = 11324
success5 = 768
success = 13580(262602)
edge swapping done!
time = 1.69058s
#v = 70106
#t = 396749
max_energy = 103.261
avg_energy = 4.00082
vertex smoothing...
success = 41952(66195)
vertex smoothing done!
time = 1.92417s
#v = 70106
#t = 396749
max_energy = 103.261
avg_energy = 3.92705
updating sclaing field ...
filter_energy = 8
is_hit_min_edge_length = 0
//////////////// pass 6 ////////////////
edge splitting...
fixed 0 tangled element
success = 31591(31591)
edge splitting done!
time = 0.355516s
#v = 101697
#t = 561544
max_energy = 101.714
avg_energy = 4.09835
edge collapsing...
fixed 0 tangled element
success(env) = 2747
success = 21128(625991)
success(env) = 181
success = 1504(208893)
success(env) = 22
success = 185(30361)
success(env) = 2
success = 36(4114)
success(env) = 0
success = 7(941)
success(env) = 0
success = 3(314)
success(env) = 0
success = 0(106)
edge collapsing done!
time = 4.80196s
#v = 78834
#t = 446410
max_energy = 31.9032
avg_energy = 3.93208
edge swapping...
fixed 0 tangled element
success3 = 1365
success4 = 11084
success5 = 733
success = 13182(299637)
edge swapping done!
time = 2.32181s
#v = 78834
#t = 445778
max_energy = 31.9032
avg_energy = 3.89862
vertex smoothing...
success = 49387(74925)
vertex smoothing done!
time = 2.07079s
#v = 78834
#t = 445778
max_energy = 24.7653
avg_energy = 3.8171
//////////////// pass 7 ////////////////
edge splitting...
fixed 0 tangled element
success = 10810(10810)
edge splitting done!
time = 0.294237s
#v = 89644
#t = 498615
max_energy = 21.6399
avg_energy = 3.87412
edge collapsing...
fixed 0 tangled element
success(env) = 1481
success = 10776(574618)
success(env) = 63
success = 411(136543)
success(env) = 6
success = 36(9337)
success(env) = 1
success = 5(864)
success(env) = 1
success = 3(138)
success(env) = 1
success = 1(92)
success(env) = 0
success = 0(16)
edge collapsing done!
time = 3.75575s
#v = 78412
#t = 443899
max_energy = 11.5316
avg_energy = 3.8237
edge swapping...
fixed 0 tangled element
success3 = 495
success4 = 6028
success5 = 306
success = 6829(285955)
edge swapping done!
time = 1.69854s
#v = 78412
#t = 443710
max_energy = 11.5316
avg_energy = 3.81349
vertex smoothing...
success = 46844(74493)
vertex smoothing done!
time = 1.97307s
#v = 78412
#t = 443710
max_energy = 11.5316
avg_energy = 3.77692
//////////////// pass 8 ////////////////
edge splitting...
fixed 0 tangled element
success = 8043(8043)
edge splitting done!
time = 0.285223s
#v = 86455
#t = 482272
max_energy = 13.0083
avg_energy = 3.822
edge collapsing...
fixed 0 tangled element
success(env) = 1166
success = 8045(559826)
success(env) = 31
success = 206(103084)
success(env) = 4
success = 21(4604)
success(env) = 0
success = 2(469)
success(env) = 0
success = 0(31)
edge collapsing done!
time = 3.37234s
#v = 78181
#t = 442590
max_energy = 11.5316
avg_energy = 3.78241
edge swapping...
fixed 0 tangled element
success3 = 230
success4 = 3810
success5 = 127
success = 4167(280351)
edge swapping done!
time = 1.61347s
#v = 78181
#t = 442487
max_energy = 11.5316
avg_energy = 3.77773
vertex smoothing...
success = 44840(74265)
vertex smoothing done!
time = 1.89036s
#v = 78181
#t = 442487
max_energy = 11.5316
avg_energy = 3.75826
updating sclaing field ...
filter_energy = 8
is_hit_min_edge_length = 0
//////////////// pass 9 ////////////////
edge splitting...
fixed 0 tangled element
success = 1219(1219)
edge splitting done!
time = 0.254546s
#v = 79400
#t = 448526
max_energy = 15.1382
avg_energy = 3.76664
edge collapsing...
fixed 0 tangled element
success(env) = 254
success = 1399(831396)
success(env) = 19
success = 68(49592)
success(env) = 2
success = 9(3217)
success(env) = 0
success = 1(741)
success(env) = 0
success = 0(118)
edge collapsing done!
time = 4.13693s
#v = 77923
#t = 441145
max_energy = 8.61096
avg_energy = 3.7607
edge swapping...
fixed 0 tangled element
success3 = 94
success4 = 1871
success5 = 78
success = 2043(275733)
edge swapping done!
time = 1.50393s
#v = 77923
#t = 441129
max_energy = 8.61096
avg_energy = 3.75924
vertex smoothing...
success = 43008(74009)
vertex smoothing done!
time = 1.83619s
#v = 77923
#t = 441129
max_energy = 8.31747
avg_energy = 3.75009
//////////////// postprocessing ////////////////
edge collapsing...
fixed 0 tangled element
success(env) = 62
success = 238(889098)
success(env) = 4
success = 31(24346)
success(env) = 1
success = 4(3308)
success(env) = 0
success = 0(386)
edge collapsing done!
time = 4.40762s
#v = 77650
#t = 439617
max_energy = 8.31747
avg_energy = 3.75069
edge collapsing...
fixed 0 tangled element
success(env) = 6760
success = 30711(573011)
success(env) = 396
success = 829(419873)
success(env) = 30
success = 64(85289)
success(env) = 3
success = 5(7999)
success(env) = 1
success = 3(575)
success(env) = 0
success = 1(284)
success(env) = 0
success = 0(106)
edge collapsing done!
time = 7.59578s
#v = 46037
#t = 258790
max_energy = 9.99983
avg_energy = 4.81014
edge swapping...
fixed 0 tangled element
success3 = 5184
success4 = 19686
success5 = 2257
success = 27127(189294)
edge swapping done!
time = 1.22345s
#v = 46037
#t = 255863
max_energy = 9.99913
avg_energy = 4.55982
edge collapsing...
fixed 0 tangled element
success(env) = 789
success = 2425(407672)
success(env) = 87
success = 140(193467)
success(env) = 1
success = 2(17341)
success(env) = 0
success = 0(342)
edge collapsing done!
time = 3.6595s
#v = 43470
#t = 240459
max_energy = 9.9993
avg_energy = 4.73164
edge swapping...
fixed 0 tangled element
success3 = 1167
success4 = 4793
success5 = 724
success = 6684(135508)
edge swapping done!
time = 0.804138s
#v = 43470
#t = 240016
max_energy = 9.99913
avg_energy = 4.67236
edge collapsing...
fixed 0 tangled element
success(env) = 89
success = 367(383709)
success(env) = 13
success = 18(42905)
success(env) = 1
success = 2(2189)
success(env) = 1
success = 1(270)
success(env) = 0
success = 0(92)
edge collapsing done!
time = 2.59632s
#v = 43082
#t = 237637
max_energy = 9.99975
avg_energy = 4.70403
edge swapping...
fixed 0 tangled element
success3 = 227
success4 = 868
success5 = 134
success = 1229(122100)
edge swapping done!
time = 0.702031s
#v = 43082
#t = 237544
max_energy = 9.99972
avg_energy = 4.69286
edge collapsing...
fixed 0 tangled element
success(env) = 9
success = 41(379983)
success(env) = 0
success = 1(5035)
success(env) = 0
success = 1(137)
success(env) = 0
success = 1(117)
success(env) = 0
success = 0(146)
edge collapsing done!
time = 2.39271s
#v = 43038
#t = 237277
max_energy = 9.99972
avg_energy = 4.69647
edge swapping...
fixed 0 tangled element
success3 = 33
success4 = 147
success5 = 22
success = 202(119555)
edge swapping done!
time = 0.668924s
#v = 43038
#t = 237266
max_energy = 9.99972
avg_energy = 4.69461
[2025-05-05 13:12:12.526] [float-tetwild] [info] mesh optimization 118.425s
[2025-05-05 13:12:12.526] [float-tetwild] [info] 
[2025-05-05 13:12:12.758] [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-05-05 13:12:14.103] [float-tetwild] [info] after winding number
[2025-05-05 13:12:14.103] [float-tetwild] [info] #v = 43038
[2025-05-05 13:12:14.106] [float-tetwild] [info] #t = 186088
[2025-05-05 13:12:14.106] [float-tetwild] [info] winding number 1.74645e+09s
[2025-05-05 13:12:14.106] [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 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

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