Meshing a packed-bed pore space with a volume-controlled Voronoi grid

Seed the interstitial space of a peclet.dem sphere packing, tessellate + clip against the sphere walls, and relax the seeds so the cell volumes follow a target — uniform, or graded to refine at the walls (an inflation layer).

voro

dem

meshing

porous-media

sdf

Author

Peclet

Published

July 4, 2026

Runs on a free Colab CPU runtime — the first cell installs peclet from PyPI.

Note

Experimental API. The SDF-walled pore-mesh functions used here — peclet.voro.sdf_voronoi_section (and optimize_pore_mesh) — are new and still evolving (see the open issues at the end). The packing (peclet.dem) is stock.

What you’ll learn

How to build a body-fitted, graded Voronoi mesh of a pore space — the Voronoi analogue of a graded CFD mesh with an inflation layer. Drop seed points into the fluid between packed spheres, tessellate, and clip every cell against the sphere surfaces (their signed-distance field) so the cells tile only the fluid. The key idea: the cell volumes are set by where the seeds go, so a good seeding heuristic — placing seeds at the target local size $s(\phi)=\mathrm{clip}(\phi, s_{lo}, s_{hi})$, small near the walls and growing into the pores — gives a graded, wall-hugging mesh directly, no volume relaxation required. You’ll compare it against a uniform mesh, cut clean cross-sections with a robust cell-section primitive, and see where a curved-wall SDF clip still leaves a hair of white.

The setup

The packing is periodic; the fluid is everything outside the spheres. The solid is described by the union signed-distance field $\phi(\mathbf{x}) = \min_i(\lVert \mathbf{x}-\mathbf{c}_i\rVert - r_i)$ ($\phi<0$ inside a sphere, $>0$ in the fluid). We seed $N$ points where $\phi>0$, build their Voronoi cells with peclet.voro, and clip each against the $\phi=0$ surfaces so the cells tile only the fluid.

A cell’s volume is roughly the cube of the local seed spacing, so to grade the mesh we place the seeds at a target size that depends on the distance to the nearest wall:

\[ s(\phi) = \mathrm{clip}(\phi,\, s_{lo},\, s_{hi}), \qquad V_{\mathrm{ref}} = s(\phi)^3 \tag{1}\]

— cells shrink to $s_{lo}$ at the walls and grow to $s_{hi}$ in the open pores. The heuristic below puts one seed roughly every $s(\phi)$, in shells around each sphere, so the mesh is the target: no optimisation needed. (If you later want the volumes made exactly equal or exactly on-target, optimize_pore_mesh relaxes the seed positions — with the caveat in the open issues.)

import importlib.util, os, subprocess, sys
_local = os.environ.get("PECLET_LOCAL_BUILD")
if _local:
    sys.path[:0] = _local.split(os.pathsep)                     # local dem/voro build(s)
elif importlib.util.find_spec("peclet") is None:
    subprocess.run([sys.executable, "-m", "pip", "install", "-q", "peclet"], check=True)

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.collections import PolyCollection
from matplotlib.patches import Circle
from matplotlib.colors import LogNorm, TwoSlopeNorm
from peclet import voro
from peclet_examples import pore_mesh as pm

plt.rcParams.update({"figure.dpi": 130, "font.size": 10, "figure.facecolor": "white",
                     "savefig.bbox": "tight"})

Kokkos::OpenMP::initialize WARNING: OMP_PROC_BIND environment variable not set
  In general, for best performance with OpenMP 4.0 or better set OMP_PROC_BIND=spread and OMP_PLACES=threads
  For best performance with OpenMP 3.1 set OMP_PROC_BIND=true
  For unit testing set OMP_PROC_BIND=false

Step 1 — pack the spheres

peclet.dem grows a random close packing (porosity ≈ 0.37).

centres, radii, L = pm.pack_spheres(n=180)
print(f"{len(centres)} spheres, box L={L:.2f}, porosity ≈ {1 - (4/3*np.pi*(radii**3).sum())/L**3:.2f}")

180 spheres, box L=5.31, porosity ≈ 0.37

Step 2 — the seeding is the volume control

The cell volumes are set by where you put the seeds, so the whole game is a good seeding heuristic. seed_graded realises a target size field $s(\phi)=\mathrm{clip}(\phi, s_{lo}, s_{hi})$ directly: it lays concentric shells around every sphere at growing distance $d$, with the radial step and the in-surface spacing both equal to the local size $s(d)$. A seed at wall-distance $\phi$ then gets a cell of size $\approx s(\phi)$ — small ($s_{lo}$) hugging the curved walls, growing to $s_{hi}$ in the open pores. That gives a body-fitted, graded, inflation-layer mesh from the seeding alone — no relaxation.

We compare two seedings: a uniform one ($s_{lo}=s_{hi}$, which meshes the walls at the bulk resolution and so can’t hug them) and a distance-graded one.

z0 = L / 2
s_uniform = pm.seed_graded(centres, radii, L, s_lo=0.19, s_hi=0.19, seed=1)   # one size everywhere
s_graded  = pm.seed_graded(centres, radii, L, s_lo=0.10, s_hi=0.36, seed=1)   # small at walls → large in pores
stages = [s_uniform, s_graded]
titles = ["uniform  (s = 0.19)", "distance-graded  (s = 0.10 → 0.36)"]
sec = [pm.section_cells(s, centres, radii, L, z0) for s in stages]   # each -> (polygons, volume, seed)
print("cells:", {t.split()[0]: len(v) for t, (_, v, _) in zip(titles, sec)})

cells: {'uniform': 557, 'distance-graded': 1535}

Step 3 — cross-section: uniform vs graded

We cut the 3-D cells with the plane $z=L/2$ and colour them by cell volume (grey discs = spheres cut by the plane). sdf_voronoi_section cuts each cell directly from its dual structure — a robust convex section that tiles the plane exactly — so no VTK and no fragile face slicing.

sl = [pm.tile_periodic(polys, vol, L) for polys, vol, _ in sec]   # +periodic images at the box faces
gv = np.asarray(sec[1][1])                                        # scale colours to the graded range
norm = LogNorm(*np.percentile(gv[gv > 0], [2, 98]))              # so its wall→pore gradient is visible

def draw(ax, polys, vals, nrm, cmap, title, win=None, lw=0.25):
    ax.add_collection(PolyCollection(polys, array=np.asarray(vals), cmap=cmap, norm=nrm,
                                     edgecolors="k", linewidths=lw, alpha=0.92))
    for c, r in zip(centres, radii):
        if abs(z0 - c[2]) < r:
            for sx in (0, -L, L):
                for sy in (0, -L, L):
                    ax.add_patch(Circle((c[0]+sx, c[1]+sy), np.sqrt(r*r-(z0-c[2])**2),
                                        facecolor="0.55", edgecolor="0.25", lw=0.5, zorder=3))
    ax.set(xlim=win[:2] if win else (0, L), ylim=win[2:] if win else (0, L), aspect="equal")
    ax.set_title(title, fontsize=10); ax.set_xticks([]); ax.set_yticks([])

fig, axes = plt.subplots(1, 2, figsize=(12, 6))
for ax, (polys, vals), t in zip(axes, sl, titles):
    draw(ax, polys, vals, norm, "viridis", t)
fig.colorbar(axes[0].collections[0], ax=axes, fraction=0.046, pad=0.02, label="cell volume (log)")
plt.show()

Figure 1: Cross-section of the two seedings, cells coloured by volume (log, shared scale). The graded mesh (right) hugs every wall with small cells and grows into the pores; the uniform mesh (left) carries one size everywhere and leaves a rim where its bulk-sized cells can’t follow the curved walls.

Both meshes fill the pore space, but only the graded one resolves the walls the way you’d want: a ring of small (dark) cells hugging every sphere, growing smoothly to large (yellow) cells in the open pores. The uniform mesh carries the bulk cell size right up to the walls, so its cells are too big to hug the curved surface — that’s the residual white rim around the spheres. The grading is not a post-process; it is exactly what the distance-shell seeding places.

Step 4 — what happens at a wall

Zooming onto one sphere makes it explicit.

cross = [(i, r*r-(z0-c[2])**2) for i, (c, r) in enumerate(zip(centres, radii)) if abs(z0-c[2]) < r]
i0 = max(cross, key=lambda t: t[1])[0]; cx, cy, w = centres[i0, 0], centres[i0, 1], 1.6
win = (cx-w, cx+w, cy-w, cy+w)
fig, ax = plt.subplots(1, 2, figsize=(11, 5.5))
for a, k, t in [(ax[0], 0, "uniform"), (ax[1], 1, "distance-graded")]:
    draw(a, *sl[k], norm, "viridis", t, win=win, lw=0.5)
plt.show()

Figure 2: Zoom at a sphere. Left: a uniform mesh — bulk-sized cells recede from the curved wall (rim). Right: the graded mesh packs small cells against the wall and grows outward — a body-fitted inflation layer.

Results — the seeding hits the target size directly

Because the shells are placed at the target size $s(\phi)$, the graded mesh’s cell volumes already match $V_{\mathrm{ref}}=s(\phi)^3$ — no relaxation needed. Colouring by $V/V_{\mathrm{ref}}$ (renormalised so 1 = on target) shows this: the graded mesh is broadly on target (pale), while the same metric against a graded target exposes how far a uniform mesh is from it.

fig, ax = plt.subplots(1, 2, figsize=(11, 5.5))
rn = TwoSlopeNorm(vmin=1/4, vcenter=1.0, vmax=4.0)
for a, s, (polys, vol, seed), t in [(ax[0], s_uniform, sec[0], "uniform mesh"),
                                    (ax[1], s_graded, sec[1], "graded mesh")]:
    vref = pm.vref_graded(pm.union_sdf(np.ascontiguousarray(s), centres, radii, L))[seed]
    relv = vol / (vref * vol.sum() / vref.sum())      # renormalise so Σ V_ref = Σ V (on-target = 1)
    draw(a, *pm.tile_periodic(polys, relv, L), rn, "coolwarm", t)
fig.colorbar(ax[1].collections[0], ax=ax, fraction=0.025, pad=0.02, label="V / V_ref (target = φ³)")
plt.show()

Figure 3: V/V_ref against the graded target φ³ (blue = smaller than target, white = on target, red = larger). The graded seeding (right) is broadly on target everywhere; the uniform mesh (left) is far off — too coarse at the walls, too fine in the pores.

Adapt this yourself

Change the grading. seed_graded(..., s_lo, s_hi) sets the wall and bulk sizes; edit the size field in pm.seed_graded from clip(φ, s_lo, s_hi) to e.g. clip(φ², …) for a thicker layer, or to a solution field (fine where a gradient is large) for solution-adaptive meshing.
Denser or coarser. Lower s_lo for a finer wall layer — cost grows like $1/s_{lo}^2$, so keep $s_{lo}\gtrsim0.08$ on a packed bed. sdf_voronoi_section cuts the section in O(N), so even a fine mesh renders quickly.
Equalise further. If you want the volumes exactly equal rather than graded, feed the seeding to voro.optimize_pore_mesh(positions, vref, …, method="steepest"|"graphamg", free_energy=…) — but see the open issue below.

Open issues (logged in ISSUES.md)

The relaxation collapses cells. optimize_pore_mesh moves seeds toward a target volume, but a position-only step can’t move seeds between pores, so from a mismatched start it collapses cells in the throats instead of redistributing. Seeding at the target size (this example) sidesteps it; a hard log-barrier keeps cells open once the start is feasible.
Curved-wall clip. The SDF clip approximates each sphere by flat (tangent) facets, so the very tightest throats keep a hair of white the mesh can’t reach; it shrinks with a finer wall layer. An exact tessellator-side curved-wall clip would close it.

(The cross-section is cut cell-by-cell with sdf_voronoi_section, which tiles exactly — the residual white is real geometry, item 2, not a slicing artefact.)

Reproduce this

pip install -e .          # + a recent peclet (peclet-voro with optimize_pore_mesh)
quarto render examples/pore-mesh-voronoi/index.qmd --execute
# local suite build instead of PyPI (dem + voro build dirs, os.pathsep-joined):
PECLET_LOCAL_BUILD=/path/to/suite/dem/build:/path/to/suite/voro/build_py \
  quarto render examples/pore-mesh-voronoi/index.qmd --execute

--- title: "Meshing a packed-bed pore space with a volume-controlled Voronoi grid" subtitle: "Seed the interstitial space of a peclet.dem sphere packing, tessellate + clip against the sphere walls, and relax the seeds so the cell volumes follow a target — uniform, or graded to refine at the walls (an inflation layer)." author: "Peclet" date: "2026-07-04" categories: [voro, dem, meshing, porous-media, sdf] jupyter: python3 --- [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/computational-chemical-engineering/peclet-examples/blob/main/examples/pore-mesh-voronoi/index.ipynb){target="_blank"}  Runs on a free Colab CPU runtime — the first cell installs `peclet` from PyPI. ::: {.callout-note} **Experimental API.** The SDF-walled pore-mesh functions used here — `peclet.voro.sdf_voronoi_section` (and `optimize_pore_mesh`) — are new and still evolving (see the open issues at the end). The packing (`peclet.dem`) is stock. ::: ## What you'll learn How to build a **body-fitted, graded Voronoi mesh of a pore space** — the Voronoi analogue of a graded CFD mesh with an inflation layer. Drop seed points into the fluid between packed spheres, tessellate, and **clip every cell against the sphere surfaces** (their signed-distance field) so the cells tile only the fluid. The key idea: the cell volumes are set by *where the seeds go*, so a good **seeding heuristic** — placing seeds at the target local size $s(\phi)=\mathrm{clip}(\phi, s_{lo}, s_{hi})$, small near the walls and growing into the pores — gives a graded, wall-hugging mesh **directly**, no volume relaxation required. You'll compare it against a uniform mesh, cut clean cross-sections with a robust cell-section primitive, and see where a curved-wall SDF clip still leaves a hair of white. ## The setup The packing is periodic; the fluid is everything outside the spheres. The solid is described by the union signed-distance field $\phi(\mathbf{x}) = \min_i(\lVert \mathbf{x}-\mathbf{c}_i\rVert - r_i)$ ($\phi<0$ inside a sphere, $>0$ in the fluid). We seed $N$ points where $\phi>0$, build their Voronoi cells with `peclet.voro`, and clip each against the $\phi=0$ surfaces so the cells tile only the fluid. A cell's volume is roughly the cube of the local seed spacing, so to grade the mesh we place the seeds at a **target size** that depends on the distance to the nearest wall: $$ s(\phi) = \mathrm{clip}(\phi,\, s_{lo},\, s_{hi}), \qquad V_{\mathrm{ref}} = s(\phi)^3 $$ {#eq-target} — cells shrink to $s_{lo}$ at the walls and grow to $s_{hi}$ in the open pores. The heuristic below puts one seed roughly every $s(\phi)$, in shells around each sphere, so the mesh *is* the target: no optimisation needed. (If you later want the volumes made *exactly* equal or exactly on-target, `optimize_pore_mesh` relaxes the seed positions — with the caveat in the open issues.) ```{python} #| label: bootstrap #| code-summary: "Environment bootstrap (installs peclet from PyPI on Colab/Binder)" import importlib.util, os, subprocess, sys _local = os.environ.get("PECLET_LOCAL_BUILD") if _local: sys.path[:0] = _local.split(os.pathsep) # local dem/voro build(s) elif importlib.util.find_spec("peclet") is None: subprocess.run([sys.executable, "-m", "pip", "install", "-q", "peclet"], check=True) import numpy as np import matplotlib.pyplot as plt from matplotlib.collections import PolyCollection from matplotlib.patches import Circle from matplotlib.colors import LogNorm, TwoSlopeNorm from peclet import voro from peclet_examples import pore_mesh as pm plt.rcParams.update({"figure.dpi": 130, "font.size": 10, "figure.facecolor": "white", "savefig.bbox": "tight"}) ``` ## Step 1 — pack the spheres `peclet.dem` grows a random close packing (porosity ≈ 0.37). ```{python} #| label: pack centres, radii, L = pm.pack_spheres(n=180) print(f"{len(centres)} spheres, box L={L:.2f}, porosity ≈ {1 - (4/3*np.pi*(radii**3).sum())/L**3:.2f}") ``` ## Step 2 — the seeding *is* the volume control The cell volumes are set by **where you put the seeds**, so the whole game is a good seeding heuristic. `seed_graded` realises a target size field $s(\phi)=\mathrm{clip}(\phi, s_{lo}, s_{hi})$ directly: it lays concentric shells around every sphere at growing distance $d$, with the radial step *and* the in-surface spacing both equal to the local size $s(d)$. A seed at wall-distance $\phi$ then gets a cell of size $\approx s(\phi)$ — small ($s_{lo}$) hugging the curved walls, growing to $s_{hi}$ in the open pores. That gives a body-fitted, graded, inflation-layer mesh **from the seeding alone — no relaxation**. We compare two seedings: a **uniform** one ($s_{lo}=s_{hi}$, which meshes the walls at the bulk resolution and so can't hug them) and a **distance-graded** one. ```{python} #| label: run z0 = L / 2 s_uniform = pm.seed_graded(centres, radii, L, s_lo=0.19, s_hi=0.19, seed=1) # one size everywhere s_graded = pm.seed_graded(centres, radii, L, s_lo=0.10, s_hi=0.36, seed=1) # small at walls → large in pores stages = [s_uniform, s_graded] titles = ["uniform (s = 0.19)", "distance-graded (s = 0.10 → 0.36)"] sec = [pm.section_cells(s, centres, radii, L, z0) for s in stages] # each -> (polygons, volume, seed) print("cells:", {t.split()[0]: len(v) for t, (_, v, _) in zip(titles, sec)}) ``` ## Step 3 — cross-section: uniform vs graded We cut the 3-D cells with the plane $z=L/2$ and colour them by cell volume (grey discs = spheres cut by the plane). `sdf_voronoi_section` cuts each cell directly from its dual structure — a robust convex section that tiles the plane exactly — so no VTK and no fragile face slicing. ```{python} #| label: fig-volume #| fig-cap: "Cross-section of the two seedings, cells coloured by volume (log, shared scale). The graded mesh (right) hugs every wall with small cells and grows into the pores; the uniform mesh (left) carries one size everywhere and leaves a rim where its bulk-sized cells can't follow the curved walls." #| fig-width: 11 sl = [pm.tile_periodic(polys, vol, L) for polys, vol, _ in sec] # +periodic images at the box faces gv = np.asarray(sec[1][1]) # scale colours to the graded range norm = LogNorm(*np.percentile(gv[gv > 0], [2, 98])) # so its wall→pore gradient is visible def draw(ax, polys, vals, nrm, cmap, title, win=None, lw=0.25): ax.add_collection(PolyCollection(polys, array=np.asarray(vals), cmap=cmap, norm=nrm, edgecolors="k", linewidths=lw, alpha=0.92)) for c, r in zip(centres, radii): if abs(z0 - c[2]) < r: for sx in (0, -L, L): for sy in (0, -L, L): ax.add_patch(Circle((c[0]+sx, c[1]+sy), np.sqrt(r*r-(z0-c[2])**2), facecolor="0.55", edgecolor="0.25", lw=0.5, zorder=3)) ax.set(xlim=win[:2] if win else (0, L), ylim=win[2:] if win else (0, L), aspect="equal") ax.set_title(title, fontsize=10); ax.set_xticks([]); ax.set_yticks([]) fig, axes = plt.subplots(1, 2, figsize=(12, 6)) for ax, (polys, vals), t in zip(axes, sl, titles): draw(ax, polys, vals, norm, "viridis", t) fig.colorbar(axes[0].collections[0], ax=axes, fraction=0.046, pad=0.02, label="cell volume (log)") plt.show() ``` Both meshes fill the pore space, but only the **graded** one resolves the walls the way you'd want: a ring of small (dark) cells hugging every sphere, growing smoothly to large (yellow) cells in the open pores. The **uniform** mesh carries the bulk cell size right up to the walls, so its cells are too big to hug the curved surface — that's the residual white rim around the spheres. The grading is not a post-process; it is exactly what the distance-shell seeding places. ## Step 4 — what happens at a wall Zooming onto one sphere makes it explicit. ```{python} #| label: fig-wall #| fig-cap: "Zoom at a sphere. Left: a uniform mesh — bulk-sized cells recede from the curved wall (rim). Right: the graded mesh packs small cells against the wall and grows outward — a body-fitted inflation layer." #| fig-width: 11 cross = [(i, r*r-(z0-c[2])**2) for i, (c, r) in enumerate(zip(centres, radii)) if abs(z0-c[2]) < r] i0 = max(cross, key=lambda t: t[1])[0]; cx, cy, w = centres[i0, 0], centres[i0, 1], 1.6 win = (cx-w, cx+w, cy-w, cy+w) fig, ax = plt.subplots(1, 2, figsize=(11, 5.5)) for a, k, t in [(ax[0], 0, "uniform"), (ax[1], 1, "distance-graded")]: draw(a, *sl[k], norm, "viridis", t, win=win, lw=0.5) plt.show() ``` ## Results — the seeding hits the target size directly Because the shells are placed *at* the target size $s(\phi)$, the graded mesh's cell volumes already match $V_{\mathrm{ref}}=s(\phi)^3$ — no relaxation needed. Colouring by $V/V_{\mathrm{ref}}$ (renormalised so 1 = on target) shows this: the graded mesh is broadly on target (pale), while the same metric against a graded target exposes how far a uniform mesh is from it. ```{python} #| label: fig-rel #| fig-cap: "V/V_ref against the graded target φ³ (blue = smaller than target, white = on target, red = larger). The graded seeding (right) is broadly on target everywhere; the uniform mesh (left) is far off — too coarse at the walls, too fine in the pores." #| fig-width: 11 fig, ax = plt.subplots(1, 2, figsize=(11, 5.5)) rn = TwoSlopeNorm(vmin=1/4, vcenter=1.0, vmax=4.0) for a, s, (polys, vol, seed), t in [(ax[0], s_uniform, sec[0], "uniform mesh"), (ax[1], s_graded, sec[1], "graded mesh")]: vref = pm.vref_graded(pm.union_sdf(np.ascontiguousarray(s), centres, radii, L))[seed] relv = vol / (vref * vol.sum() / vref.sum()) # renormalise so Σ V_ref = Σ V (on-target = 1) draw(a, *pm.tile_periodic(polys, relv, L), rn, "coolwarm", t) fig.colorbar(ax[1].collections[0], ax=ax, fraction=0.025, pad=0.02, label="V / V_ref (target = φ³)") plt.show() ``` ## Adapt this yourself - **Change the grading.** `seed_graded(..., s_lo, s_hi)` sets the wall and bulk sizes; edit the size field in `pm.seed_graded` from `clip(φ, s_lo, s_hi)` to e.g. `clip(φ², …)` for a thicker layer, or to a *solution* field (fine where a gradient is large) for solution-adaptive meshing. - **Denser or coarser.** Lower `s_lo` for a finer wall layer — cost grows like $1/s_{lo}^2$, so keep $s_{lo}\gtrsim0.08$ on a packed bed. `sdf_voronoi_section` cuts the section in O(N), so even a fine mesh renders quickly. - **Equalise further.** If you want the volumes *exactly* equal rather than graded, feed the seeding to `voro.optimize_pore_mesh(positions, vref, …, method="steepest"|"graphamg", free_energy=…)` — but see the open issue below. ## Open issues (logged in ISSUES.md) 1. **The relaxation collapses cells.** `optimize_pore_mesh` moves seeds toward a target volume, but a position-only step can't move seeds *between* pores, so from a mismatched start it collapses cells in the throats instead of redistributing. Seeding *at* the target size (this example) sidesteps it; a hard log-barrier keeps cells open once the start is feasible. 2. **Curved-wall clip.** The SDF clip approximates each sphere by flat (tangent) facets, so the very tightest throats keep a hair of white the mesh can't reach; it shrinks with a finer wall layer. An exact tessellator-side curved-wall clip would close it. *(The cross-section is cut cell-by-cell with `sdf_voronoi_section`, which tiles exactly — the residual white is real geometry, item 2, not a slicing artefact.)* ## Reproduce this ```bash pip install -e . # + a recent peclet (peclet-voro with optimize_pore_mesh) quarto render examples/pore-mesh-voronoi/index.qmd --execute # local suite build instead of PyPI (dem + voro build dirs, os.pathsep-joined): PECLET_LOCAL_BUILD=/path/to/suite/dem/build:/path/to/suite/voro/build_py \ quarto render examples/pore-mesh-voronoi/index.qmd --execute ```