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