8 · Contact: a Ball Pushing into a Soft Block#
Contact is just another energy: a penalty spring that switches on wherever a vertex pokes inside an obstacle. We use a ball and the signed distance
(negative = inside = penalized).
We do four things:
a static ball probe — the block stays frozen while a ghost ball slides through it, and we read off the contact energy vs the ball’s position,
a static ground probe — a rigid ball descends into a flat ground and we read off the plane-contact energy, which peaks when the ball is halfway through,
a dynamic simulation — a strong ball that visibly dents a hanging, pinned block,
a stiffer-contact rerun — the same scripted ball, with a 10× firmer penalty.
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import simkit
from simkit.solvers import newton_solver
import utils
The block, the ball, and the contact energy#
utils.SphereContact(K, r, M_n, dim) wraps the penalty that sums over every
vertex inside the ball; utils.PlaneContact(K, point, normal, M_n, dim) does the
same against a flat half-space. We build the hanging block and pin its top edge.
X, T = utils.triangulated_grid(nx=10, ny=12, width=0.6, height=0.8)
X[:, 1] += 0.4 # hang it: y in [0, 0.8]
n, dim = X.shape
pin_idx = np.where(X[:, 1] >= X[:, 1].max() - 1e-6)[0] # top edge pinned
BALL_R = 0.16
K_CONTACT = 5e4 # strong contact
The block simulator#
BlockSim is the dynamic simulator for the soft block. It composes the standard
term helpers from utils — the Neo-Hookean material, a pin spring,
gravity, the sphere contact penalty, and Backward-Euler inertia —
and sums them one per line in energy / gradient / hessian. set_ball(c)
re-aims the contact obstacle, and step() advances one Backward-Euler step with
Newton. The same class drives both the soft-contact run (Experiment 3) and the
stiffer-contact run (Experiment 4); only K_contact changes.
class BlockSim:
"""Dynamic soft block under gravity, pinned at top, denting against a ball."""
def __init__(self, X, T, K_contact, ball_r, h=0.02):
self.p = utils.precompute(X, T, ym=200.0, pr=0.4, rho=1.0, gravity=-4.0)
self.psi = utils.make_material("Neo-Hookean")
self.pin = utils.PenaltySpring(self.p.n, self.p.dim, 1e6).set(pin_idx, X[pin_idx])
self.gravity = utils.Gravity(self.p.f_g)
self.contact = utils.SphereContact(K_contact, ball_r, self.p.M_n, self.p.dim)
self.inertia = utils.Inertia(self.p.M, h)
self.x = X.reshape(-1, 1).astype(float)
self.v = np.zeros_like(self.x)
def set_ball(self, c):
self.contact.set_center(c)
return self
def energy(self, x):
E_elastic = self.psi.energy(x, self.p)
E_pin = self.pin.energy(x)
E_gravity = self.gravity.energy(x)
E_contact = self.contact.energy(x)
E_inertia = self.inertia.energy(x)
return E_elastic + E_pin + E_gravity + E_contact + E_inertia
def gradient(self, x):
g_elastic = self.psi.gradient(x, self.p)
g_pin = self.pin.gradient(x)
g_gravity = self.gravity.gradient(x)
g_contact = self.contact.gradient(x)
g_inertia = self.inertia.gradient(x)
return g_elastic + g_pin + g_gravity + g_contact + g_inertia
def hessian(self, x):
H_elastic = self.psi.hessian(x, self.p)
H_pin = self.pin.hessian(x)
H_contact = self.contact.hessian(x)
H_inertia = self.inertia.hessian(x)
return H_elastic + H_pin + H_contact + H_inertia
def step(self):
self.inertia.update(self.x, self.v)
x_new = newton_solver(self.x, self.energy, self.gradient, self.hessian,
max_iter=100, do_line_search=True)
self.v = (x_new - self.x) / self.inertia.h
self.x = x_new
return self.x.reshape(self.p.n, self.p.dim)
Experiment 1 — static probe: the block does not move#
Freeze the block at its rest pose and slide a ghost ball up through it. We evaluate the contact energy of that static block as a function of the ball’s position — no stepping, no deformation. The signed-distance background shows the contact region (blue); the energy on the right rises exactly when the ball overlaps more of the block.
ball_y = np.linspace(-0.5, 0.55, 44)
centers = [np.array([0.0, cy]) for cy in ball_y]
X_flat = X.reshape(-1, 1)
# block is STATIC: evaluate contact energy of the rest block for each ball position
probe_energy = np.array([utils.SphereContact(K_CONTACT, BALL_R, simkit.massmatrix(X, T, rho=1.0), dim)
.set_center(c).energy(X_flat) for c in centers])
static_block = [X.copy() for _ in centers] # block never deforms here
fig, anim = utils.animate_dynamics(static_block, T, ball_y, {"contact": probe_energy},
lims=((-0.6, 0.6), (-0.7, 0.95)), xlabel="ball height", ylabel="contact energy",
colors=utils.ENERGY_COLORS, ball_centers=centers, ball_radius=BALL_R, sdf=True,
scene_title="STATIC block, ghost ball sweeping (blue = inside)",
title="contact energy vs ball position", fps=18)
utils.show_video(fig, anim, "media/08_ball_probe.mp4", fps=18)
Experiment 2 — ground probe: energy peaks halfway through#
Same idea, but now the obstacle is the ground. We drop a rigid ball toward a flat plane and read off the plane-contact energy of the ball’s vertices. As the ball sinks in, more of it goes below the surface and the energy climbs — reaching its maximum when the ball is halfway through (its center on the ground), then falling as it rises back out. The blue background is the ground’s signed distance
(negative = inside the ground). The ball here is not simulated; it is a probe.
ball_X, ball_T = utils.ball_mesh_2d(radius=0.22, n_segments=48)
ball_M = simkit.massmatrix(ball_X, ball_T, rho=1.0)
GROUND_Y = 0.0
floor_p, floor_n = np.array([0.0, GROUND_Y]), np.array([0.0, 1.0])
K_GROUND, R2 = 5e3, 0.22
ground = utils.PlaneContact(K_GROUND, floor_p, floor_n, ball_M, dim)
# sweep the ball center from above the ground down to ground level (halfway
# submerged) and back up; energy peaks at the deepest point
cy = np.concatenate([np.linspace(R2 + 0.2, GROUND_Y, 22), np.linspace(GROUND_Y, R2 + 0.2, 22)])
ball_states = [ball_X + np.array([0.0, y]) for y in cy]
ground_energy = np.array([ground.energy(s.reshape(-1, 1)) for s in ball_states])
print(f"peak ground-contact energy at ball-center y = {cy[int(np.argmax(ground_energy))]:.2f} (ground at {GROUND_Y})")
fig, anim = utils.animate_dynamics(ball_states, ball_T, np.arange(len(ball_states)),
{"contact": ground_energy}, lims=((-0.6, 0.6), (-0.45, 0.7)),
xlabel="frame", ylabel="ground-contact energy", colors=utils.ENERGY_COLORS,
floor_y=GROUND_Y, floor_sdf=True, mesh_face=utils.BALL_FACE, mesh_edge=utils.BALL_EDGE,
scene_title="rigid ball sinking into the ground (blue = inside)",
title="energy peaks at halfway through", fps=18)
utils.show_video(fig, anim, "media/08_ground_probe.mp4", fps=18)
peak ground-contact energy at ball-center y = 0.00 (ground at 0.0)
Experiment 3 — dynamic simulation: the ball dents the block#
Now let the block respond. Each ball position is a Backward-Euler step that minimizes elastic + pin + gravity + strong contact + inertia. The block visibly deforms as the ball drives up into it.
sim = BlockSim(X, T, K_CONTACT, BALL_R, h=0.02)
ball_y = np.concatenate([np.linspace(-0.5, 0.4, 22), np.linspace(0.4, -0.5, 22)])
centers, states, contact_E = [], [], []
for cy in ball_y:
c_ball = np.array([0.0, cy])
sim.set_ball(c_ball)
U = sim.step()
centers.append(c_ball)
states.append(U.copy())
contact_E.append(sim.contact.energy(sim.x))
print(f"max block displacement: {max(np.abs(s - X).max() for s in states):.3f}")
fig, anim = utils.animate_dynamics(states, T, np.arange(len(states)), {"contact": np.array(contact_E)},
lims=((-0.6, 0.6), (-0.7, 0.95)), xlabel="frame", ylabel="contact energy",
colors=utils.ENERGY_COLORS, ball_centers=centers, ball_radius=BALL_R, sdf=True,
scene_title="ball denting the block (simulated)", title="contact energy", fps=18)
utils.show_video(fig, anim, "media/08_contact_ball.mp4", fps=18)
max block displacement: 0.278
Experiment 4 — stiffer contact: the same ball, a firmer penalty#
We repeat Experiment 3 exactly — same block, same pin, same gravity, and
the same scripted ball motion — but crank the contact stiffness \(k\)
up by 10× (K_CONTACT_STIFF = 5e5). We do not rewrite any physics: we just
build a second BlockSim with the higher K_contact. A stiffer penalty
spring fights penetration harder: vertices are allowed to poke far less into the
ball, so the block is shoved away more decisively. Watch the dent grow and the
contact-energy curve climb well above the soft-contact run.
K_CONTACT_STIFF = 5e5 # 10x stiffer contact than Experiment 3
sim_stiff = BlockSim(X, T, K_CONTACT_STIFF, BALL_R, h=0.02)
ball_y = np.concatenate([np.linspace(-0.5, 0.4, 22), np.linspace(0.4, -0.5, 22)]) # SAME scripted motion
centers, states_stiff, contact_E_stiff = [], [], []
for cy in ball_y:
c_ball = np.array([0.0, cy])
sim_stiff.set_ball(c_ball)
U = sim_stiff.step()
centers.append(c_ball)
states_stiff.append(U.copy())
contact_E_stiff.append(sim_stiff.contact.energy(sim_stiff.x))
print(f"max block displacement (stiff contact): {max(np.abs(s - X).max() for s in states_stiff):.3f}")
fig, anim = utils.animate_dynamics(states_stiff, T, np.arange(len(states_stiff)), {"contact": np.array(contact_E_stiff)},
lims=((-0.6, 0.6), (-0.7, 0.95)), xlabel="frame", ylabel="contact energy",
colors=utils.ENERGY_COLORS, ball_centers=centers, ball_radius=BALL_R, sdf=True,
scene_title="stiffer contact (10x): ball denting the block", title="contact energy (stiff)", fps=18)
utils.show_video(fig, anim, "media/08_contact_ball_stiff.mp4", fps=18)
max block displacement (stiff contact): 0.490
Takeaways#
Penalty contact is an energy that turns on where the signed distance is negative; minimizing it pushes penetrating vertices back out.
The static probe shows the energy purely as a function of overlap (no simulation); the dynamic run lets the block deform in response.
Deeper penetration → quadratically larger contact energy.
Contact stiffness is a knob: a stiffer penalty (Experiment 4) tolerates less penetration, so the same scripted ball dents the block more and the contact energy spikes higher — at the cost of a stiffer, harder-to-solve system.