Custom-shaped particles: pack a boxful of beans (peclet.dem)

Author a particle as a signed distance function, turn it into a DEM particle — surface shell, moment of inertia and all — and pack a boxful.

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

What you’ll learn

peclet.dem is not limited to spheres. A particle is a signed distance field (SDF) plus a set of surface points, and the collision solver works for any shape you can describe implicitly. In this example you will:

1. Author a particle — a kidney bean — as a signed distance function.
2. Turn it into a DEM particle with peclet.dem.build_particle, which samples the SDF onto a grid, extracts a surface point shell by marching cubes, and computes the mass, centre of mass and full inertia tensor by voxel integration — returning the particle in its principal-axis frame so the solver’s inertia is exact (no hand-tuned numbers).
3. Pack a boxful of beans by Lubachevsky–Stillinger growth and measure the packing fraction.

Everything runs on a CPU in well under a minute.

import importlib.util, os, subprocess, sys
_local = os.environ.get("PECLET_LOCAL_BUILD")
if _local:
    for p in _local.split(os.pathsep):
        sys.path.insert(0, p)
elif importlib.util.find_spec("peclet") is None:
    subprocess.run([sys.executable, "-m", "pip", "install", "-q", "peclet"], check=True)
# build_particle turns an SDF into a surface shell via marching cubes (scikit-image);
# the packing is shown as an interactive 3-D figure with plotly.
for _mod, _pip in (("skimage", "scikit-image"), ("plotly", "plotly")):
    if importlib.util.find_spec(_mod) is None:
        subprocess.run([sys.executable, "-m", "pip", "install", "-q", _pip], check=True)

import numpy as np
import math, time
import matplotlib.pyplot as plt
from peclet import dem
from peclet.dem import build_particle

plt.rcParams.update({"figure.dpi": 130, "font.size": 10, "axes.axisbelow": True,
                     "figure.facecolor": "white", "savefig.bbox": "tight"})
print("peclet.dem backend:", dem.execution_space)

peclet.dem backend: Cuda

Step 1 — Author a bean as an SDF

An SDF is a function that returns, for any point, the signed distance to the surface — negative inside the solid, positive outside. A kidney bean is a bent, tapering, slightly flattened sausage: sweep a sphere along a circular arc, shrink its radius toward the ends, and squash it a little in z. Each of those is one line acting on an (N, 3) array of points:

def bean(arc_r=0.55, arc_a=1.15, r0=0.30, taper=0.55, zflat=0.72):
    """Kidney bean: a sphere of radius r0 swept along an arc of radius arc_r spanning
    +/- arc_a radians, tapering toward the tips and flattened in z."""
    def f(p):
        x, y, z = p[:, 0], p[:, 1], p[:, 2]
        theta = np.clip(np.arctan2(y, x), -arc_a, arc_a)   # angle along the spine
        cx, cy = arc_r * np.cos(theta), arc_r * np.sin(theta)   # nearest point on the arc
        dist = np.sqrt((x - cx)**2 + (y - cy)**2 + (z / zflat)**2)
        radius = r0 * (1.0 - taper * (np.abs(theta) / arc_a)**2)  # taper toward the tips
        return dist - radius
    return f

Tip

Any callable f(points) -> distances works — a bent capsule like this, boolean CSG of primitives, or the output of a dedicated implicit-modelling library such as fogleman/sdf. We keep it to a few NumPy lines here so the example is self-contained.

Step 2 — Build the DEM particle

build_particle does the heavy lifting: it samples the SDF on a grid, extracts a surface point shell (marching cubes, thinned to a controlled density), and integrates the mass properties. It returns the shape in its principal-axis body frame, so the diagonal inverse inertia the solver uses is exact.

shape = build_particle(
    bean(),
    bounds=((0.0, -0.7, -0.4), (0.95, 0.7, 0.4)),   # a box enclosing the solid
    resolution=72,
    target_shell_points=220,
)
print(f"mass            = {shape.mass:.4f}")
print(f"bounding radius = {shape.bounding_radius:.3f}")
print(f"surface points  = {len(shape.shell)}")
print(f"principal inertia (unit mass) = {np.round(shape.inertia / shape.mass, 4)}")
print(f"  -> anisotropy I_max / I_min = {shape.inertia.max() / shape.inertia.min():.2f}"
      "  (a bean tumbles differently about each axis)")

mass            = 0.1863
bounding radius = 0.687
surface points  = 300
principal inertia (unit mass) = [0.0342 0.0991 0.1155]
  -> anisotropy I_max / I_min = 3.37  (a bean tumbles differently about each axis)

The particle is fully described by its SDF (the field, whose zero level set is the surface) and the surface shell (the collision probes marched off it). The xy slice shows the kidney profile; the shell is the point cloud the solver tests:

fig = plt.figure(figsize=(9, 3.8))

ax = fig.add_subplot(121)
mid = shape.grid.shape[2] // 2
o, sp, g = shape.origin, shape.spacing, shape.grid
extent = [o[0], o[0] + sp[0]*(g.shape[0]-1), o[1], o[1] + sp[1]*(g.shape[1]-1)]
im = ax.imshow(g[:, :, mid].T, origin="lower", extent=extent, cmap="RdBu", vmin=-0.35, vmax=0.35)
ax.contour(g[:, :, mid].T, levels=[0], colors="k", linewidths=1.3, extent=extent)
ax.set_title("SDF slice (z = 0)"); ax.set_aspect("equal")
fig.colorbar(im, ax=ax, shrink=0.75, label="signed distance")

ax = fig.add_subplot(122, projection="3d")
sh = shape.shell
ax.scatter(sh[:, 0], sh[:, 1], sh[:, 2], s=5, c=sh[:, 1], cmap="viridis")
ax.set_title(f"surface point shell ({len(sh)} pts)"); ax.set_box_aspect((1, 1.3, 0.8))
fig.tight_layout(); plt.show()

Step 3 — Pack a boxful

We pack the beans by the Lubachevsky–Stillinger protocol: the particles start small (so there is no initial overlap), then grow while a thermostat keeps the assembly fluid. Growth is gated on the measured overlap, so they settle into gentle contact rather than being forced to interpenetrate. A final quench freezes the packing.

The particle’s shape, surface shell and inertia are shared by every particle via shape.apply_to(sim); per-particle position, orientation and scale are yours to set — exactly as for the built-in sphere.

N = 80
phi_target = 0.45
D = (N * shape.mass / phi_target) ** (1/3)      # cube sized so full-size beans give phi_target
plen = 2 * shape.bounding_radius

sim = dem.Simulation(N)
shape.apply_to(sim)                              # <-- the imported bean becomes THE shape
sim.set_domain((0, 0, 0), (D, D, D))
sim.enable_periodicity(False, False, False)
for point, normal in [((0,0,0),(1,0,0)), ((D,0,0),(-1,0,0)), ((0,0,0),(0,1,0)),
                      ((0,D,0),(0,-1,0)), ((0,0,0),(0,0,1)), ((0,0,D),(0,0,-1))]:
    sim.add_plane(point, normal)                 # six walls
sim.set_gravity(0, 0, 0)
sim.set_material_params(0.3, 0.3, 0.0)
sim.set_solver_iterations(40, 8)

rng = np.random.default_rng(11)
p = rng.uniform(0.15*D, 0.85*D, (N, 4)).astype(np.float32); p[:, 3] = 1.0
sim.set_positions(p)
sim.set_velocities(rng.normal(0, 0.4, (N, 3)).astype(np.float32))
q = rng.normal(0, 1, (N, 4)).astype(np.float32); q /= np.linalg.norm(q, axis=1, keepdims=True)
sim.set_quaternions(q)                           # random initial orientations
sim.set_angular_velocities(np.zeros((N, 3), np.float32))

rate, dt, crit = 1.0, 0.005, 0.02 * plen
sim.set_scales(np.full(N, 1.0, np.float32))
sim.set_growth_params(rate, 0.06)
sim.set_thermostat(0.6, 0.004)
t0 = time.time()
for i in range(1200):                            # grow to full size, gated on overlap
    grow = sim.get_max_overlap() < crit and float(sim.get_scales().mean()) < 0.999
    sim.set_growth_params(rate if grow else 0.0, sim.get_growth_factor())
    sim.step(dt)
sim.set_thermostat(0, 1); sim.set_material_params(0, 0, 0)
for i in range(300):                             # quench to a static packing
    sim.step(dt)

phi = phi_target * float(np.mean(sim.get_scales() ** 3))
print(f"packed {N} beans in {time.time()-t0:.1f} s")
print(f"packing fraction phi ~ {phi:.3f}")
print(f"max overlap = {sim.get_max_overlap()/plen*100:.1f}% of the bean length")

packed 80 beans in 3.8 s
packing fraction phi ~ 0.450
max overlap = 0.0% of the bean length

The packing — as solid beans

For a proper picture we don’t need the point shell. shape.grid is the bean’s SDF in its body frame, so its zero level set is the surface: one marching-cubes call gives a triangle mesh — the same shape.vertices / shape.faces you would write to an STL — that we place, rotate and scale for every bean. Rendered with Plotly the figure is interactive: drag to rotate, scroll to zoom.

import plotly.graph_objects as go
from skimage.measure import marching_cubes

pos = sim.get_positions().reshape(-1, 3)
quat = sim.get_quaternions().reshape(-1, 4)          # (x, y, z, w)
scale = sim.get_scales().ravel()

# bean surface mesh in its body frame (coarsened for a light interactive figure; == shape.vertices/faces -> STL)
stride = 5
grid = np.ascontiguousarray(shape.grid[::stride, ::stride, ::stride])
cverts, cfaces, _, _ = marching_cubes(grid, level=0.0, spacing=tuple(np.array(shape.spacing) * stride))
cverts = cverts + shape.origin

def qrot(q, V):                                      # rotate verts by quaternion q = (x, y, z, w)
    qv = q[:3]; t = 2.0 * np.cross(np.broadcast_to(qv, V.shape), V)
    return V + q[3] * t + np.cross(np.broadcast_to(qv, V.shape), t)

V, F, C, off = [], [], [], 0
for i in range(N):
    w = qrot(quat[i], cverts * scale[i]) + pos[i]    # place, rotate and scale this bean
    V.append(w); F.append(cfaces + off); C.append(np.full(len(w), pos[i, 2])); off += len(w)
V, F, C = np.concatenate(V), np.concatenate(F), np.concatenate(C)

fig = go.Figure(go.Mesh3d(
    x=V[:, 0], y=V[:, 1], z=V[:, 2], i=F[:, 0], j=F[:, 1], k=F[:, 2],
    intensity=C, colorscale="Viridis", showscale=False, flatshading=True,
    lighting=dict(ambient=0.5, diffuse=0.8, specular=0.2, roughness=0.6),
    lightposition=dict(x=D, y=1.6 * D, z=1.3 * D)))
fig.update_layout(width=560, height=560, margin=dict(l=0, r=0, t=34, b=0),
                  title=f"N = {N} beans, φ ≈ {phi:.2f} — drag to rotate",
                  scene=dict(aspectmode="data", xaxis=dict(range=[0, D]), yaxis=dict(range=[0, D]),
                             zaxis=dict(range=[0, D]), xaxis_title="x", yaxis_title="y", zaxis_title="z"))
fig.show()

(a) The packed beans in 3-D (drag to rotate). Every bean is its own mesh — marched off the shape’s body-frame SDF (a ready-made STL) — placed at the position, orientation and grown size the DEM solver assigns it, coloured by height.

(b)

Figure 1

Adapt this yourself

• Change the particle. Swap bean() for any SDF. Boolean CSG works too — a dumbbell is a smooth union of two lobes and a neck:

  def sphere(c, r): c = np.asarray(c); return lambda p: np.linalg.norm(p - c, axis=1) - r
  def smooth_union(a, b, k=0.08):
      def f(p):
          da, db = a(p), b(p); h = np.clip(0.5 + 0.5*(db - da)/k, 0, 1)
          return db*(1-h) + da*h - k*h*(1-h)
      return f
  dumbbell = smooth_union(sphere((-0.4,0,0), 0.32), sphere((0.4,0,0), 0.32))
  shape = build_particle(dumbbell, ((-0.8,-0.4,-0.4), (0.8,0.4,0.4)), resolution=64)

  Convex shapes (rounded cubes, ellipsoids) resolve contacts most cleanly; gentle concavities like the bean’s inner curve are fine. For a full CSG algebra, feed build_particle a callable from fogleman/sdf.

• Real mass / density. apply_to uses unit mass by default. For a real density, call shape.apply_to(sim, unit_mass=False) and set sim.set_inv_mass(1/mass) — the builder exposes shape.mass and shape.inertia at your chosen density=.

• Export the particle / packing. shape.to_stl("bean.stl") writes the bean’s surface mesh to STL (also available as shape.vertices / shape.faces) — open it in Blender/ParaView/Ovito, or render the packing offline for a publication figure. sim.write_vtp("beans.vtp") writes particle centres/orientations for ParaView/Ovito; sim.get_sdf_grid((128,128,128)) reconstructs the pore-space SDF of the whole bed — the bridge to a CFD run with peclet.flow (see the random packed bed example).

• Go bigger / on a GPU. Increase N; on a CUDA/HIP build the same script runs on the device. A distributed build adds sim.init_mpi(...) / sim.step_mpi(...).

Reproduce this

# from a checkout of peclet-examples
pip install -e '.[sim]'          # or: pip install peclet scikit-image
# refresh the rendered page from a local build of the suite:
PECLET_LOCAL_BUILD=/path/to/suite/dem/build_omp \
  quarto render examples/sdf-particle-packing/index.qmd --execute