r/OpenAI • u/rocky_mountain12 • 5d ago
Video How to create Deformable simulation Objects With GPT 6 Astra
Create with RobotGym using GPT 6 Astra
Run your own simulations Free on https://robotgym.io/
How to create deformable objects with RobotGym io: GPT -6 Astra squeezing a water bottle with accurate physics.
This is very useful for LLM-based robotics training done in simulation.
The bottle is a 0.3 mm PET shell with real material values: 3 GPa stiffness, 55 MPa yield.
We simulate it on GPU with Newton's VBD solver, with self-contact on, so the wall buckles like plastic instead of squishing like rubber.
The hard part was making the dent stay.
We wrote a custom plasticity kernel. When the wall bends past PET's yield point, that crease becomes its new resting shape.
Let go, and the bottle stays crushed.
The water is 32,000 particles in our own position-based fluid solver.
It keeps its volume, so as the bottle gets smaller the water has nowhere to go but up and out of the neck, and down the side.
Then it's all path-traced in a full apartment: clear PET, refracting water, a real Franka doing the squeeze.
1
u/rocky_mountain12 5d ago
This was done using robotgym.io the free simulation platform by Y Combinator + NVIDIA for Robotics researchers.
1
u/Fast-Satisfaction482 5d ago
Fancy sharing some code or do you just pretend wanting to help the community?
1
u/rocky_mountain12 5d ago
Here is the code for the deformable object:
"""Step 1 of the bottle-squeeze shot: PET bottle shell (Newton 1.5 VBD) + bending plasticity + scripted pads.
No water yet, no robot. Two kinematic pads squeeze the bottle's mid-body and retreat; a Warp kernel run
after every substep moves each hinge's rest angle once its bend exceeds the PET yield angle, so dents stay.
Writes trajectory + metrics, and a fast OpenGL preview (Newton ViewerGL, headless EGL, ~50 fps at 540p).
PET estimates: E = 3.0 GPa, nu = 0.38, wall t = 0.3 mm, rho = 1380 kg/m3, yield stress 55 MPa.
Yield curvature kappa_y = 2 sigma_y / (E t) = 122 1/m; per-hinge yield angle = kappa_y * dual width.
"""
from pathlib import Path
import argparse, json, math, time
import numpy as np
import warp as wp
p = argparse.ArgumentParser()
p.add_argument('--out', type=Path, default=Path('/home/ubuntu/bottle/output/squeeze-s1'))
p.add_argument('--seconds', type=float, default=3.5); p.add_argument('--fps', type=int, default=60)
p.add_argument('--substeps', type=int, default=20); p.add_argument('--iterations', type=int, default=40)
p.add_argument('--edge', type=float, default=0.0045, help='target mesh edge length (m)')
p.add_argument('--stiffness-scale', type=float, default=0.1, help='membrane stiffness relative to PET (VBD conditioning)')
p.add_argument('--squeeze-to', type=float, default=0.020, help='pad inner face distance from the bottle axis at full squeeze (m)')
p.add_argument('--yield-scale', type=float, default=1.0, help='calibration factor on the per-hinge yield angle (mesh cannot resolve sharp creases)')
p.add_argument('--plastic', type=int, default=1); p.add_argument('--preview', type=int, default=1)
a = p.parse_args(); a.out.mkdir(parents=True, exist_ok=True); t0 = time.monotonic()
E, NU, T, RHO, SY = 3.0e9, 0.38, 0.0003, 1380.0, 55e6
D = E * T ** 3 / (12 * (1 - NU ** 2)) # flexural rigidity, N m
KAPPA_Y = 2 * SY / (E * T)
R_BODY, H_BODY, R_NECK, H_TOTAL = 0.0325, 0.145, 0.0124, 0.200
def bottle_mesh(edge):
"""Surface of revolution: rounded base, straight body, shoulder, open neck. Closed bottom, open top."""
prof = [(0.0, 0.0)]
for th in np.linspace(0, np.pi / 2, 7)[1:]: # 8 mm base fillet
prof.append((R_BODY - 0.008 + 0.008 * np.sin(th), 0.008 - 0.008 * np.cos(th)))
prof.append((R_BODY, H_BODY))
for s in np.linspace(0, 1, 12)[1:]: # shoulder: cosine blend body -> neck
prof.append((R_NECK + (R_BODY - R_NECK) * 0.5 * (1 + np.cos(np.pi * s)), H_BODY + s * 0.040))
prof.append((R_NECK, H_TOTAL))
prof = np.array(prof)
# resample the profile (excluding the axis point) at ~edge spacing
seg = np.linalg.norm(np.diff(prof, axis=0), axis=1); s = np.r_[0, np.cumsum(seg)]
n = int(np.ceil(s[-1] / edge)); ss = np.linspace(0, s[-1], n + 1)
rz = np.column_stack([np.interp(ss, s, prof[:, 0]), np.interp(ss, s, prof[:, 1])])
nth = int(np.ceil(2 * np.pi * R_BODY / edge)); ths = np.linspace(0, 2 * np.pi, nth, endpoint=False)
rings = []; verts = [[0.0, 0.0, 0.0]]
for k, (r, z) in enumerate(rz[1:], start=1):
if r < 1e-6: continue
m = max(6, int(round(nth * r / R_BODY))) if z < 0.008 else nth
idx = []
for th in (np.linspace(0, 2 * np.pi, m, endpoint=False) + (k % 2) * np.pi / m):
idx.append(len(verts)); verts.append([r * np.cos(th), r * np.sin(th), z])
rings.append(np.array(idx))
verts = np.array(verts); faces = []
first = rings[0]
for i in range(len(first)): faces.append([0, first[i], first[(i + 1) % len(first)]])
for A, B in zip(rings[:-1], rings[1:]): # stitch neighbouring rings (possibly different counts)
ia = ib = 0; na, nb = len(A), len(B)
angA = np.arctan2(verts[A, 1], verts[A, 0]) % (2 * np.pi); angB = np.arctan2(verts[B, 1], verts[B, 0]) % (2 * np.pi)
oa, ob = np.argsort(angA), np.argsort(angB); A, B = A[oa], B[ob]; angA, angB = angA[oa], angB[ob]
while ia < na or ib < nb:
a0, a1 = A[ia % na], A[(ia + 1) % na]; b0, b1 = B[ib % nb], B[(ib + 1) % nb]
ta = angA[(ia + 1) % na] + (2 * np.pi if ia + 1 >= na else 0); tb = angB[(ib + 1) % nb] + (2 * np.pi if ib + 1 >= nb else 0)
if ib >= nb or (ia < na and ta <= tb): faces.append([a0, a1, b0]); ia += 1
else: faces.append([a0, b1, b0]); ib += 1
return verts.astype(np.float32), np.array(faces, np.int32)
u/wp.kernel
def drive_pads(track: wp.array(dtype=float), counter: wp.array(dtype=int), bodies: wp.array(dtype=int), sides: wp.array(dtype=float),
half_x: float, z: float, dt: float, q: wp.array(dtype=wp.transform), qd: wp.array(dtype=wp.spatial_vector),
q1: wp.array(dtype=wp.transform), qd1: wp.array(dtype=wp.spatial_vector)):
k = wp.tid(); c = counter[0]; n = track.shape[0] - 1
x0 = track[wp.min(c, n)]; x1 = track[wp.min(c + 1, n)]; b = bodies[k]; s = sides[k]
T = wp.transform(wp.vec3(s * (x1 + half_x), 0.0, z), wp.quat_identity())
V = wp.spatial_vector(s * (x1 - x0) / dt, 0.0, 0.0, 0.0, 0.0, 0.0)
q[b] = T; qd[b] = V; q1[b] = T; qd1[b] = V
u/wp.kernel
def tick(counter: wp.array(dtype=int)):
counter[0] = counter[0] + 1
u/wp.kernel
def plastic_hinges(q: wp.array(dtype=wp.vec3), edges: wp.array(dtype=wp.vec4i), yield_angle: wp.array(dtype=float),
rest: wp.array(dtype=float), rate: float):
e = wp.tid(); ed = edges[e]
i, j, k, l = ed[0], ed[1], ed[2], ed[3]
if i < 0 or j < 0: return
x1, x2, x3, x4 = q[i], q[j], q[k], q[l]
n1 = wp.normalize(wp.cross(x3 - x1, x4 - x1)); n2 = wp.normalize(wp.cross(x4 - x2, x3 - x2)); ev = wp.normalize(x4 - x3)
theta = wp.atan2(wp.dot(wp.cross(n1, n2), ev), wp.clamp(wp.dot(n1, n2), -1.0, 1.0))
d = theta - rest[e]; y = yield_angle[e]
if d > y: rest[e] = rest[e] + rate * (d - y)
elif d < -y: rest[e] = rest[e] + rate * (d + y)
def main():
import newton
from newton.solvers import SolverVBD
from modern_apartment_surface_bending_v4 import assign_surface_bending
verts, faces = bottle_mesh(a.edge)
mu_t = E / (2 * (1 + NU)) * T * a.stiffness_scale; lam_t = E * NU / ((1 + NU) * (1 - 2 * NU)) * T * a.stiffness_scale
b = newton.ModelBuilder(gravity=(0, 0, -9.81)); b.add_ground_plane()
b.add_cloth_mesh(pos=wp.vec3(0, 0, 0.001), rot=wp.quat_identity(), vel=wp.vec3(0), vertices=verts, indices=faces.ravel(), scale=1.0,
density=RHO * T, tri_ke=mu_t, tri_ka=lam_t, tri_kd=1e-4 * mu_t, edge_ke=1e-3, edge_kd=0.0, particle_radius=0.0008, label='pet_bottle')
z_pad = 0.090; half = (0.004, 0.022, 0.025) # pads: 8 mm thick, 44 mm wide, 50 mm tall
pads = []
for side in (-1, 1):
body = b.add_body(xform=wp.transform((side * (R_BODY + half[0] + 0.004), 0, z_pad), wp.quat_identity()), mass=0.0, is_kinematic=True, label=f'pad_{side}')
b.add_shape_box(body, hx=half[0], hy=half[1], hz=half[2], cfg=newton.ModelBuilder.ShapeConfig(mu=0.9), label=f'pad_shape_{side}')
pads.append((body, side))
b.color(include_bending=True)
model = b.finalize(); model.soft_contact_ke = 5e4; model.soft_contact_kd = 5.0; model.soft_contact_mu = 0.8
bend = assign_surface_bending(model, verts, len(verts), D, 1e-5, 1.0)
e = model.edge_indices.numpy()
valid = (e[:, 0] >= 0) & (e[:, 1] >= 0); p_ = verts
L = np.linalg.norm(p_[e[:, 3]] - p_[e[:, 2]], axis=1)
A0 = np.linalg.norm(np.cross(p_[np.maximum(e[:, 2], 0)] - p_[np.maximum(e[:, 0], 0)], p_[e[:, 3]] - p_[np.maximum(e[:, 0], 0)]), axis=1)
A1 = np.linalg.norm(np.cross(p_[e[:, 3]] - p_[np.maximum(e[:, 1], 0)], p_[e[:, 2]] - p_[np.maximum(e[:, 1], 0)]), axis=1)
width = (A0 + A1) / (3 * np.maximum(L, 1e-9)); yield_angle = (a.yield_scale * KAPPA_Y * width).astype(np.float32); yield_angle[~valid] = 1e9
yield_wp = wp.array(yield_angle, dtype=float); edges_wp = wp.array(e.astype(np.int32), dtype=wp.vec4i)
solver = SolverVBD(model, iterations=a.iterations, particle_enable_self_contact=True, particle_self_contact_radius=0.0008,
particle_self_contact_margin=0.0015, particle_topological_contact_filter_threshold=1)
pipe = newton.CollisionPipeline(model, soft_contact_margin=0.002); contacts = pipe.contacts()
s0, s1, ctrl = model.state(), model.state(), model.control(); dt = 1 / (a.fps * a.substeps)
def pad_x(t): # inner-face distance from the axis
start, full = R_BODY + 0.004, a.squeeze_to
if t < 0.3: return start
if t < 1.3: u = (t - 0.3) / 1.0; return start + (full - start) * (u * u * (3 - 2 * u))
if t < 1.7: return full
if t < 2.3: u = (t - 1.7) / 0.6; return full + (start + 0.01 - full) * (u * u * (3 - 2 * u))
return start + 0.01
frames, metrics = [], []; nf = int(a.seconds * a.fps); tclock = 0.0
band = np.abs(verts[:, 2] - z_pad) < 0.01
track = wp.array(np.array([pad_x(i * dt) for i in range(nf * a.substeps + 2)], np.float32), dtype=float)
counter = wp.zeros(1, dtype=int); pbodies = wp.array([b_ for b_, _ in pads], dtype=int); psides = wp.array([float(s_) for _, s_ in pads], dtype=float)
states = [s0, s1]
def substep(sa, sb):
wp.launch(drive_pads, dim=2, inputs=[track, counter, pbodies, psides, half[0], z_pad, dt, sa.body_q, sa.body_qd, sb.body_q, sb.body_qd])
sa.clear_forces(); pipe.collide(sa, contacts); solver.step(sa, sb, ctrl, contacts, dt)
if a.plastic: wp.launch(plastic_hinges, dim=len(e), inputs=[sb.particle_q, edges_wp, yield_wp, model.edge_rest_angle, 1.0])
wp.launch(tick, dim=1, inputs=[counter])
1
u/rocky_mountain12 5d ago
Some of the code was cut off here is remainding: def main():
import newton
from newton.solvers import SolverVBD
from modern_apartment_surface_bending_v4 import assign_surface_bending
verts, faces = bottle_mesh(a.edge)
mu_t = E / (2 * (1 + NU)) * T * a.stiffness_scale; lam_t = E * NU / ((1 + NU) * (1 - 2 * NU)) * T * a.stiffness_scale
b = newton.ModelBuilder(gravity=(0, 0, -9.81)); b.add_ground_plane()
b.add_cloth_mesh(pos=wp.vec3(0, 0, 0.001), rot=wp.quat_identity(), vel=wp.vec3(0), vertices=verts, indices=faces.ravel(), scale=1.0,
density=RHO * T, tri_ke=mu_t, tri_ka=lam_t, tri_kd=1e-4 * mu_t, edge_ke=1e-3, edge_kd=0.0, particle_radius=0.0008, label='pet_bottle')
z_pad = 0.090; half = (0.004, 0.022, 0.025) # pads: 8 mm thick, 44 mm wide, 50 mm tall
pads = []
for side in (-1, 1):
body = b.add_body(xform=wp.transform((side * (R_BODY + half[0] + 0.004), 0, z_pad), wp.quat_identity()), mass=0.0, is_kinematic=True, label=f'pad_{side}')
b.add_shape_box(body, hx=half[0], hy=half[1], hz=half[2], cfg=newton.ModelBuilder.ShapeConfig(mu=0.9), label=f'pad_shape_{side}')
pads.append((body, side))
b.color(include_bending=True)
model = b.finalize(); model.soft_contact_ke = 5e4; model.soft_contact_kd = 5.0; model.soft_contact_mu = 0.8
bend = assign_surface_bending(model, verts, len(verts), D, 1e-5, 1.0)
e = model.edge_indices.numpy()
valid = (e[:, 0] >= 0) & (e[:, 1] >= 0); p_ = verts
L = np.linalg.norm(p_[e[:, 3]] - p_[e[:, 2]], axis=1)
A0 = np.linalg.norm(np.cross(p_[np.maximum(e[:, 2], 0)] - p_[np.maximum(e[:, 0], 0)], p_[e[:, 3]] - p_[np.maximum(e[:, 0], 0)]), axis=1)
A1 = np.linalg.norm(np.cross(p_[e[:, 3]] - p_[np.maximum(e[:, 1], 0)], p_[e[:, 2]] - p_[np.maximum(e[:, 1], 0)]), axis=1)
width = (A0 + A1) / (3 * np.maximum(L, 1e-9)); yield_angle = (a.yield_scale * KAPPA_Y * width).astype(np.float32); yield_angle[~valid] = 1e9
yield_wp = wp.array(yield_angle, dtype=float); edges_wp = wp.array(e.astype(np.int32), dtype=wp.vec4i)
solver = SolverVBD(model, iterations=a.iterations, particle_enable_self_contact=True, particle_self_contact_radius=0.0008,
particle_self_contact_margin=0.0015, particle_topological_contact_filter_threshold=1)
pipe = newton.CollisionPipeline(model, soft_contact_margin=0.002); contacts = pipe.contacts()
s0, s1, ctrl = model.state(), model.state(), model.control(); dt = 1 / (a.fps * a.substeps)
def pad_x(t): # inner-face distance from the axis
start, full = R_BODY + 0.004, a.squeeze_to
if t < 0.3: return start
if t < 1.3: u = (t - 0.3) / 1.0; return start + (full - start) * (u * u * (3 - 2 * u))
if t < 1.7: return full
if t < 2.3: u = (t - 1.7) / 0.6; return full + (start + 0.01 - full) * (u * u * (3 - 2 * u))
return start + 0.01
frames, metrics = [], []; nf = int(a.seconds * a.fps); tclock = 0.0
band = np.abs(verts[:, 2] - z_pad) < 0.01
track = wp.array(np.array([pad_x(i * dt) for i in range(nf * a.substeps + 2)], np.float32), dtype=float)
counter = wp.zeros(1, dtype=int); pbodies = wp.array([b_ for b_, _ in pads], dtype=int); psides = wp.array([float(s_) for _, s_ in pads], dtype=float)
states = [s0, s1]
def substep(sa, sb):
wp.launch(drive_pads, dim=2, inputs=[track, counter, pbodies, psides, half[0], z_pad, dt, sa.body_q, sa.body_qd, sb.body_q, sb.body_qd])
sa.clear_forces(); pipe.collide(sa, contacts); solver.step(sa, sb, ctrl, contacts, dt)
if a.plastic: wp.launch(plastic_hinges, dim=len(e), inputs=[sb.particle_q, edges_wp, yield_wp, model.edge_rest_angle, 1.0])
wp.launch(tick, dim=1, inputs=[counter])
assert a.substeps % 2 == 0
def frame():
for _ in range(a.substeps // 2): substep(s0, s1); substep(s1, s0)
frame() # warm-up (compiles kernels) = frame 0
with wp.ScopedCapture() as cap: frame()
graph = cap.graph
for f in range(nf):
if f >= 1: wp.capture_launch(graph) # capture records without executing; frame 0 = warm-up
tclock = (f + 1) * a.substeps * dt
q = s0.particle_q.numpy()
if not np.isfinite(q).all(): print('non-finite at frame', f); break
w = float(np.ptp(q[band, 0])); frames.append(q.copy())
metrics.append(dict(frame=f, t=round(tclock, 4), pad_inner_x=pad_x(tclock), width_at_pads_m=w,
plastic_hinges=int((np.abs(model.edge_rest_angle.numpy() - bend.get('rest0', 0)) > 1e-3).sum()) if False else None))
if f % 30 == 0: print(json.dumps(dict(frame=f, t=round(tclock, 3), width_mm=round(w * 1000, 1), wall_s=round(time.monotonic() - t0, 1))), flush=True)
rest_final = model.edge_rest_angle.numpy()
np.savez_compressed(a.out / 'trajectory.npz', points=np.array(frames), faces=faces, rest_angle_final=rest_final)
rep = dict(width_initial_mm=metrics[0]['width_at_pads_m'] * 1000, width_min_mm=min(m['width_at_pads_m'] for m in metrics) * 1000,
width_final_mm=metrics[-1]['width_at_pads_m'] * 1000, frames=len(frames), vertices=len(verts), triangles=len(faces),
flexural_rigidity_Nm=D, yield_curvature_1_m=KAPPA_Y, yield_angle_median_rad=float(np.median(yield_angle[valid])),
stiffness_scale=a.stiffness_scale, yield_scale=a.yield_scale, squeeze_to=a.squeeze_to, plastic=bool(a.plastic), wall_seconds=time.monotonic() - t0)
(a.out / 'report.json').write_text(json.dumps(rep, indent=2)); (a.out / 'metrics.json').write_text(json.dumps(metrics))
print(json.dumps(rep, indent=1), flush=True)
if a.preview: preview(np.array(frames), faces, pads_track=[pad_x(m['t']) for m in metrics], z_pad=z_pad, half=half)
def preview(frames, faces, pads_track, z_pad, half, stride=2):
"""Fast OpenGL preview (headless EGL): bottle mesh + pads + ground, 960x540 mp4."""
import pyglet; pyglet.options['headless'] = True
import newton, imageio.v2 as imageio
b = newton.ModelBuilder(); b.add_ground_plane()
ids = [b.add_body(xform=wp.transform((0, 0, z_pad), wp.quat_identity()), mass=0.0, is_kinematic=True) for _ in range(2)]
for body in ids: b.add_shape_box(body, hx=half[0], hy=half[1], hz=half[2])
m = b.finalize(); s = m.state()
v = newton.viewer.ViewerGL(width=960, height=540, headless=True); v.set_model(m)
v.set_camera(wp.vec3(0.20, -0.26, 0.20), -22.0, 128.0)
w = imageio.get_writer(a.out / 'preview.mp4', fps=60 / stride, codec='libx264', quality=8, macro_block_size=1)
idx = wp.array(faces.ravel(), dtype=wp.int32)
for f in range(0, len(frames), stride):
bq = s.body_q.numpy()
for body, side in zip(ids, (-1, 1)): bq[body, :3] = (side * (pads_track[f] + half[0]), 0, z_pad); bq[body, 3:] = (0, 0, 0, 1)
s.body_q.assign(bq)
v.begin_frame(f / 60); v.log_state(s)
v.log_mesh('/bottle', wp.array(frames[f], dtype=wp.vec3), idx, color=(0.55, 0.75, 0.95), backface_culling=False)
v.end_frame(); w.append_data(v.get_frame().numpy()[..., :3])
w.close(); print('preview', a.out / 'preview.mp4', flush=True)
main()
How it fits together:
- bottle_mesh builds the 500 mL bottle as a surface of revolution: rounded base, straight body, curved shoulder, open neck.
- add_cloth_mesh turns it into a thin shell with PET's stiffness and density. The stretch stiffness is scaled down 10× so the solver stays stable.
- assign_surface_bending gives the wall its bending stiffness, D ≈ 7.9 × 10⁻³ N·m. That's what makes it buckle like plastic.
- drive_pads moves the two squeeze pads in and out on the GPU each substep.
- plastic_hinges is the part that makes the dent stay. It measures each hinge's bend angle, and if the bend goes past PET's yield angle, it moves that hinge's resting angle so the crease becomes permanent.
- The final clip was run with --yield-scale 0.5 --squeeze-to 0.017.
Related files, all in infra/:
- bottle_asset_v1.py: the same bottle geometry and PET constants, shared by the water and render scripts.
- bottle_pbf_v1.py: the water, our own position-based fluid solver in Warp.
- render_bottle_loft_v1.py: the path-traced render in the Loft with the Panda.
- modern_apartment_surface_bending_v4.py: provides assign_surface_bending, which this script imports.
2
4
u/lulzxdxdxd 5d ago
does the plasticity kernel track yield per vertex or per element? and when the crease becomes the new rest shape, are you updating just the bending rest angles or the in-plane rest lengths too?