{ "cells": [ { "cell_type": "markdown", "id": "cell-title", "metadata": {}, "source": [ "# Simplicits with MPM (Material Point Method)\n", "\n", "This notebook demonstrates one-way coupling between Simplicits soft-body simulation\n", "and MPM granular materials (sand). Two deformable cubes act as static obstacles inside\n", "the MPM model — each frame, Simplicits particles are mirrored into the MPM state via\n", "`move_static_particles_kernel` so the sand sees the current cube positions." ] }, { "cell_type": "markdown", "id": "cell-imports-md", "metadata": {}, "source": [ "## Imports\n", "\n", "We import Kaolin, Newton, Warp, and **k3d** for interactive 3D visualization.\n", "`SolverImplicitMPM` drives the sand simulation; `SimplicitsModelBuilder` /\n", "`SimplicitsSolver` drive the soft-body cubes. `ParticleFlags` marks particles\n", "as active or static." ] }, { "cell_type": "code", "execution_count": null, "id": "0ba3fb6d-3d19-4944-ab62-a33d9c3b43f7", "metadata": {}, "outputs": [], "source": [ "!pip install -q newton[examples]==1.2.0" ] }, { "cell_type": "code", "execution_count": null, "id": "cell-imports", "metadata": {}, "outputs": [], "source": [ "import os\n", "import sys\n", "import threading\n", "from pathlib import Path\n", "\n", "import numpy as np\n", "import torch\n", "import warp as wp\n", "from IPython.display import display\n", "from ipywidgets import Button, VBox\n", "\n", "import newton\n", "from newton import ModelBuilder, Contacts\n", "from newton.solvers import SolverImplicitMPM\n", "from newton._src.geometry import ParticleFlags\n", "\n", "import kaolin\n", "import kaolin.physics.simplicits\n", "from kaolin.experimental.newton.builder import SimplicitsModelBuilder\n", "from kaolin.experimental.newton.solver import SimplicitsSolver\n", "sys.path.append(str(Path(\"..\")))\n", "from tutorial_common import COMMON_DATA_DIR\n", "\n", "try:\n", " import k3d\n", "except ModuleNotFoundError:\n", " !pip install -q k3d\n", " import k3d\n", " raise RuntimeError(\"k3d has been installed. Do a hard refresh of the page (Ctrl + Shift + R) to avoid issue with k3d\")" ] }, { "cell_type": "markdown", "id": "cell-config-md", "metadata": {}, "source": [ "## Configuration\n", "\n", "Key constants:\n", "- **`SOFT_YOUNGS_MODULUS` / `POISSON_RATIO` / `DENSITY`** — elastic material parameters.\n", "- **`FRAME_DT`** — simulation timestep. MPM runs `sim_substeps=10` sub-steps per frame.\n", "- **`COLLISION_PARTICLE_RADIUS`** — radius used for Simplicits self-collision detection.\n", "- **`MPM_VOXEL_SIZE`** — grid cell size for the MPM solver.\n", "- **`SAND_BOUNDS_LO/HI`** — bounding box for sand particle spawn region.\n", "- **`CUBE_TRANSFORM_1/2`** — initial 4×4 transforms for the two soft cubes." ] }, { "cell_type": "code", "execution_count": null, "id": "cell-config", "metadata": {}, "outputs": [], "source": [ "# Simulation parameters\n", "FRAME_DT = 1.0 / 100.0\n", "FLOOR_PLANE = 0.0\n", "\n", "# Simplicits parameters\n", "SOFT_YOUNGS_MODULUS = 1e4\n", "POISSON_RATIO = 0.45\n", "DENSITY = 500\n", "APPROX_VOLUME = 0.5\n", "NUM_SAMPLES = 10000\n", "\n", "# Contact parameters\n", "SELF_CONTACT_RADIUS = 0.01\n", "SELF_CONTACT_MARGIN = 0.01\n", "SOFT_CONTACT_KE = 100\n", "SOFT_CONTACT_KD = 2e-3\n", "SELF_CONTACT_FRICTION = 0.25\n", "SOFT_CONTACT_MAX = 1000000\n", "\n", "# Collision parameters\n", "COLLISION_PARTICLE_RADIUS = 0.05\n", "DETECTION_RATIO = 1.5\n", "IMPENETRABLE_BARRIER_RATIO = 0.25\n", "COLLISION_PENALTY = 1000.0\n", "MAX_CONTACT_PAIRS = 10000\n", "COLLISION_FRICTION = 0.5\n", "\n", "# MPM parameters\n", "MPM_VOXEL_SIZE = 0.07\n", "MPM_TOLERANCE = 1.0e-6\n", "MPM_MAX_ITERATIONS = 250\n", "MPM_PARTICLE_KE = 1.0e15\n", "MPM_PARTICLE_KD = 0.01\n", "MPM_PARTICLE_MU = 0.5\n", "PARTICLES_PER_CELL = 3\n", "\n", "# Inter-model collision parameters\n", "MAX_INTER_CONTACTS = 20000\n", "INTER_COLLISION_STIFFNESS = 1e4\n", "INTER_COLLISION_DAMPING = 100.0\n", "CONTACT_FORCE_SCALE = 1.0\n", "MAX_VELOCITY_CHANGE = 10.0\n", "\n", "# Sand spawn parameters\n", "SAND_BOUNDS_LO = np.array([0.25, -0.5, 1.5])\n", "SAND_BOUNDS_HI = np.array([0.75, 0.5, 1.75])\n", "SAND_DENSITY = 2500.0\n", "\n", "# Simplicits spawn transforms\n", "CUBE_TRANSFORM_1 = torch.tensor([\n", " [1, 0, 0, 0.4],\n", " [0, 1, 0, 0.1],\n", " [0, 0, 1, 0.3],\n", " [0, 0, 0, 1]\n", "], dtype=torch.float32, device='cuda')\n", "\n", "CUBE_TRANSFORM_2 = torch.tensor([\n", " [1, 0, 0, 0.4],\n", " [0, 1, 0, 0.1],\n", " [0, 0, 1, 1.0],\n", " [0, 0, 0, 1]\n", "], dtype=torch.float32, device='cuda')" ] }, { "cell_type": "markdown", "id": "cell-kernels-md", "metadata": {}, "source": [ "## Warp Kernels\n", "**`move_static_particles_kernel`** — each frame, copies Simplicits particle positions\n", " and velocities into the MPM static obstacle slots." ] }, { "cell_type": "code", "execution_count": null, "id": "cell-kernels", "metadata": {}, "outputs": [], "source": [ "wp.clear_kernel_cache()\n", "\n", "@wp.kernel\n", "def move_static_particles_kernel(\n", " static_idx: wp.array(dtype=int),\n", " particle_q: wp.array(dtype=wp.vec3),\n", " particle_qd: wp.array(dtype=wp.vec3),\n", " sim_pts: wp.array(dtype=wp.vec3),\n", " sim_vel: wp.array(dtype=wp.vec3),\n", "):\n", " \"\"\"Update static MPM particles to mirror current Simplicits particle positions.\"\"\"\n", " tid = wp.tid()\n", " st_idx = static_idx[tid]\n", " particle_q[st_idx] = sim_pts[tid]\n", " particle_qd[st_idx] = sim_vel[tid]" ] }, { "cell_type": "markdown", "id": "cell-simplicits-md", "metadata": {}, "source": [ "## Simplicits Object and Model\n", "\n", "`get_simplicits_cube()` loads a unit cube mesh, scales it by 0.5, samples interior\n", "points, and creates a affinely deformable `SimplicitsObject`. It returns both `sim_obj` (for physics)\n", "and `mesh` (for k3d visualization).\n", "\n", "Two cubes are added at `CUBE_TRANSFORM_1` and `CUBE_TRANSFORM_2`. Their object IDs\n", "(`obj_id_0`, `obj_id_1`) are captured for later use in `get_object_deformed_pts`." ] }, { "cell_type": "code", "execution_count": null, "id": "cell-simplicits", "metadata": {}, "outputs": [], "source": [ "def get_simplicits_cube():\n", " \"\"\"Load cube mesh and create a single handle Simplicits object. Returns (sim_obj, mesh).\"\"\"\n", " mesh_path = os.path.join(COMMON_DATA_DIR, \"meshes\")\n", " mesh = kaolin.io.import_mesh(\n", " mesh_path + \"/cube.obj\", triangulate=True\n", " ).to('cuda')\n", " mesh.vertices = kaolin.ops.pointcloud.center_points(\n", " mesh.vertices.unsqueeze(0), normalize=True\n", " ).squeeze(0) * 0.5\n", " orig_vertices = mesh.vertices.clone()\n", "\n", " uniform_pts = torch.rand(NUM_SAMPLES, 3, device='cuda') * (\n", " orig_vertices.max(dim=0).values - orig_vertices.min(dim=0).values\n", " ) + orig_vertices.min(dim=0).values\n", "\n", " physics_points = kaolin.physics.simplicits.PhysicsPoints(\n", " pts=uniform_pts,\n", " yms=SOFT_YOUNGS_MODULUS,\n", " prs=POISSON_RATIO,\n", " rhos=DENSITY,\n", " appx_vol=APPROX_VOLUME\n", " ).cuda()\n", "\n", " sim_obj = kaolin.physics.simplicits.SimplicitsObject.create_rigid(\n", " physics_points=physics_points\n", " ) # create_rigid constructs an *deformable* object with a single handle for fast construction. At higher YM values it acts as rigid.\n", " return sim_obj, mesh\n", "\n", "\n", "# Create Simplicits objects\n", "sim_obj, mesh = get_simplicits_cube()\n", "orig_vertices = mesh.vertices.clone()\n", "\n", "# Build Simplicits model\n", "simplicits_builder = SimplicitsModelBuilder()\n", "simplicits_builder.add_ground_plane()\n", "\n", "simplicits_builder.add_simplicits_object(sim_obj, num_qp=2000, init_transform=CUBE_TRANSFORM_1, renderable_pts=orig_vertices.clone().detach())\n", "simplicits_builder.add_simplicits_object(sim_obj, num_qp=2000, init_transform=CUBE_TRANSFORM_2, renderable_pts=orig_vertices.clone().detach())\n", "\n", "simplicits_builder.add_simplicits_collisions(\n", " collision_particle_radius=COLLISION_PARTICLE_RADIUS,\n", " detection_ratio=DETECTION_RATIO,\n", " impenetrable_barrier_ratio=IMPENETRABLE_BARRIER_RATIO,\n", " collision_penalty=COLLISION_PENALTY,\n", " max_contact_pairs=MAX_CONTACT_PAIRS,\n", " friction=COLLISION_FRICTION\n", ")\n", "\n", "simplicits_model = simplicits_builder.finalize()\n", "simplicits_state_0 = simplicits_model.state()\n", "simplicits_state_1 = simplicits_model.state()\n", "\n", "simplicits_solver = SimplicitsSolver(simplicits_model)\n", "simplicits_model.simplicits_scene.max_newton_steps = 5\n", "simplicits_model.simplicits_scene.max_ls_steps = 10\n", "simplicits_model.simplicits_scene.newton_hessian_regularizer = 1e-4\n", "\n", "print(f\"Simplicits model: {simplicits_model.particle_count} particles\")" ] }, { "cell_type": "markdown", "id": "cell-mpm-md", "metadata": {}, "source": [ "## MPM Setup\n", "\n", "The MPM model holds two types of particles:\n", "\n", "1. **Sand particles** (indices `0 : mpm_static_start_idx`) — active granular material\n", " spawned in `SAND_BOUNDS_LO/HI` with `PARTICLES_PER_CELL` per voxel cell.\n", "2. **Static Simplicits obstacle particles** (indices `mpm_static_start_idx :`) — copies\n", " of Simplicits integration points added with `flags=0` (inactive for MPM dynamics).\n", " Each frame, `move_static_particles_kernel` updates their positions to match the\n", " current Simplicits state so the MPM solver sees updated obstacles." ] }, { "cell_type": "code", "execution_count": null, "id": "cell-mpm", "metadata": {}, "outputs": [], "source": [ "def _spawn_particles(builder, voxel_size, bounds_lo, bounds_hi, density, flags):\n", " \"\"\"Spawn particles on a regular grid within a bounding box.\"\"\"\n", " res = np.array(\n", " np.ceil(PARTICLES_PER_CELL * (bounds_hi - bounds_lo) / voxel_size), dtype=int\n", " )\n", " px = np.linspace(bounds_lo[0], bounds_hi[0], res[0] + 1)\n", " py = np.linspace(bounds_lo[1], bounds_hi[1], res[1] + 1)\n", " pz = np.linspace(bounds_lo[2], bounds_hi[2], res[2] + 1)\n", " points = np.stack(np.meshgrid(px, py, pz)).reshape(3, -1).T\n", " cell_size = (bounds_hi - bounds_lo) / res\n", " cell_volume = np.prod(cell_size)\n", " radius = np.max(cell_size) * 1.0\n", " print(f\"mpm particleradius: {radius}\")\n", " mass = np.prod(cell_volume) * density\n", " rng = np.random.default_rng(42)\n", " points += 2.0 * radius * (rng.random(points.shape) - 0.5)\n", " first_id = len(builder.particle_q)\n", " builder.add_particles(\n", " pos=[p for p in points],\n", " vel=[(0.0, 0.0, 0.0)] * points.shape[0],\n", " mass=[mass] * points.shape[0],\n", " radius=[radius] * points.shape[0],\n", " flags=[flags] * points.shape[0],\n", " )\n", " return np.arange(first_id, first_id + points.shape[0], dtype=int)\n", "\n", "\n", "def emit_particles(builder, voxel_size):\n", " \"\"\"Spawn sand particles in SAND_BOUNDS region.\"\"\"\n", " return _spawn_particles(\n", " builder, voxel_size, SAND_BOUNDS_LO, SAND_BOUNDS_HI, SAND_DENSITY, ParticleFlags.ACTIVE\n", " )\n", "\n", "\n", "sim_substeps = 10\n", "frame_dt = FRAME_DT\n", "sim_dt = frame_dt / sim_substeps\n", "\n", "mpm_builder = ModelBuilder()\n", "SolverImplicitMPM.register_custom_attributes(mpm_builder)\n", "mpm_builder.add_ground_plane()\n", "\n", "sand_particles_ids = emit_particles(mpm_builder, voxel_size=MPM_VOXEL_SIZE)\n", "\n", "# Mirror Simplicits particles as static obstacles in MPM\n", "sim_pts = simplicits_model.sim_z_to_full(simplicits_state_0.sim_z).numpy()\n", "mpm_static_start_idx = len(mpm_builder.particle_q)\n", "\n", "mpm_builder.add_particles(\n", " pos=[p for p in sim_pts],\n", " vel=[(0.0, 0.0, 0.0)] * sim_pts.shape[0],\n", " mass=[1.0] * sim_pts.shape[0],\n", " radius=[0.1] * sim_pts.shape[0],\n", " flags=[0] * sim_pts.shape[0],\n", ")\n", "\n", "mpm_static_idx_end = len(mpm_builder.particle_q)\n", "mpm_static_idx = wp.array(\n", " np.arange(mpm_static_start_idx, mpm_static_idx_end, dtype=int),\n", " dtype=wp.int32, device='cuda'\n", ")\n", "\n", "mpm_model = mpm_builder.finalize()\n", "\n", "mpm_options = SolverImplicitMPM.Config()\n", "mpm_options.voxel_size = MPM_VOXEL_SIZE\n", "mpm_options.tolerance = MPM_TOLERANCE\n", "mpm_options.max_iterations = MPM_MAX_ITERATIONS\n", "\n", "mpm_state_0 = mpm_model.state()\n", "mpm_state_1 = mpm_model.state()\n", "\n", "mpm_solver = SolverImplicitMPM(mpm_model, mpm_options)\n", "\n", "mpm_model.particle_ke = MPM_PARTICLE_KE\n", "mpm_model.particle_kd = MPM_PARTICLE_KD\n", "mpm_model.particle_mu = MPM_PARTICLE_MU\n", "\n", "print(f\"MPM model: {mpm_model.particle_count} particles \"\n", " f\"({mpm_static_start_idx} sand + {mpm_static_idx_end - mpm_static_start_idx} static)\")" ] }, { "cell_type": "markdown", "id": "cell-simloop-md", "metadata": {}, "source": [ "## Simulation Loop\n", "\n", "Each call to `simulate()` advances the scene by one frame:\n", "\n", "1. **Simplicits step** — one Newton solve advances the two soft cubes by `frame_dt`.\n", "2. **Static particle update** — `move_static_particles_kernel` copies the new Simplicits\n", " particle positions/velocities into the MPM static obstacle slots.\n", "3. **MPM substeps** — `sim_substeps=10` implicit MPM steps advance the sand by `sim_dt` each." ] }, { "cell_type": "code", "execution_count": null, "id": "cell-simloop", "metadata": {}, "outputs": [], "source": [ "from ipywidgets import Output\n", "out = Output()\n", "\n", "def simulate():\n", " with out:\n", " global simplicits_state_0, simplicits_state_1, mpm_state_0, mpm_state_1\n", " \n", " # Simplicits step (one big step per frame)\n", " simplicits_state_0.clear_forces()\n", " simplicits_state_1.clear_forces()\n", " contacts = simplicits_model.collide(simplicits_state_0)\n", " simplicits_solver.step(\n", " simplicits_state_0, simplicits_state_1,\n", " None, contacts, frame_dt\n", " )\n", " simplicits_state_0, simplicits_state_1 = simplicits_state_1, simplicits_state_0\n", " \n", " # Update MPM static obstacle particles to match current Simplicits positions\n", " wp_sim_pts = simplicits_model.sim_z_to_full(simplicits_state_1.sim_z)\n", " wp_sim_vel = simplicits_model.sim_z_dot_to_full(simplicits_state_0.sim_z_dot)\n", " wp.launch(\n", " kernel=move_static_particles_kernel,\n", " dim=mpm_static_idx.shape[0],\n", " inputs=[\n", " mpm_static_idx,\n", " mpm_state_0.particle_q,\n", " mpm_state_0.particle_qd,\n", " wp_sim_pts,\n", " wp_sim_vel,\n", " ],\n", " )\n", "\n", " # MPM substeps\n", " for _ in range(sim_substeps):\n", " mpm_state_0.clear_forces()\n", " mpm_state_1.clear_forces()\n", " mpm_solver.step(mpm_state_0, mpm_state_1, None, None, sim_dt)\n", " mpm_state_0, mpm_state_1 = mpm_state_1, mpm_state_0" ] }, { "cell_type": "markdown", "id": "cell-viewer-md", "metadata": {}, "source": [ "## Visualization with k3d\n", "\n", "Two k3d objects are rendered:\n", "\n", "- **`mesh_soft`** (blue) — the two Simplicits cubes combined into a single mesh.\n", " Updated via `get_object_deformed_pts` each frame.\n", "- **`points_sand`** (yellow) — the active sand particles\n", " (`mpm_state_0.particle_q[:mpm_static_start_idx]`).\n", " The static Simplicits obstacle copies are hidden.\n", "\n", "Play/Stop/Reset buttons control the background simulation thread." ] }, { "cell_type": "code", "execution_count": null, "id": "cell-viewer", "metadata": {}, "outputs": [], "source": [ "orig_vertices_np = orig_vertices.cpu().numpy().astype(np.float32)\n", "orig_faces_np = mesh.faces.cpu().numpy().astype(np.uint32)\n", "num_verts = orig_vertices_np.shape[0]\n", "\n", "combined_faces_np = np.concatenate(\n", " [orig_faces_np, orig_faces_np + num_verts]\n", ").astype(np.uint32)\n", "\n", "plot = k3d.plot(camera_auto_fit=False)\n", "plot.camera = [3, 3, 3, 0, 0, 0, 0, 0, 1]\n", "\n", "mesh_soft = k3d.mesh(\n", " np.concatenate([orig_vertices_np, orig_vertices_np]).astype(np.float32),\n", " combined_faces_np,\n", " color=0x3399ff\n", ")\n", "\n", "sand_pts_np = mpm_state_0.particle_q.numpy()[:mpm_static_start_idx].astype(np.float32)\n", "points_sand = k3d.points(sand_pts_np, color=0xffcc00, point_size=0.05)\n", "\n", "plot += mesh_soft\n", "plot += points_sand\n", "\n", "# Floor plane at FLOOR_PLANE height\n", "_s = 3.0 # half-extent\n", "floor_verts = np.array([\n", " [-_s, -_s, FLOOR_PLANE],\n", " [ _s, -_s, FLOOR_PLANE],\n", " [ _s, _s, FLOOR_PLANE],\n", " [-_s, _s, FLOOR_PLANE],\n", "], dtype=np.float32)\n", "floor_faces = np.array([[0, 1, 2], [0, 2, 3]], dtype=np.uint32)\n", "mesh_floor = k3d.mesh(floor_verts, floor_faces, color=0xcccccc, opacity=0.5)\n", "plot += mesh_floor\n", "\n", "plot.display()\n", "\n", "\n", "def update_vis_mesh():\n", " # Soft body cubes\n", " with out:\n", " print(\".\")\n", " simplicits_model.simplicits_scene.sim_z = simplicits_state_0.sim_z\n", " v0 = simplicits_model.simplicits_scene.get_object_deformed_pts(0, 'rendered')\n", " v1 = simplicits_model.simplicits_scene.get_object_deformed_pts(1, 'rendered')\n", " mesh_soft.vertices = np.concatenate(\n", " [v0.cpu().numpy(), v1.cpu().numpy()]\n", " ).astype(np.float32)\n", " # Sand particles (active only)\n", " points_sand.positions = mpm_state_0.particle_q.numpy()[:mpm_static_start_idx].astype(np.float32)\n", "\n", "\n", "sim_running = [False]\n", "sim_thread = [None]\n", "\n", "\n", "def run_sim_loop():\n", " while sim_running[0]:\n", " simulate()\n", " update_vis_mesh()\n", "\n", "\n", "def on_play(b):\n", " if not sim_running[0]:\n", " sim_running[0] = True\n", " sim_thread[0] = threading.Thread(target=run_sim_loop, daemon=True)\n", " sim_thread[0].start()\n", "\n", "\n", "def on_stop(b):\n", " sim_running[0] = False\n", "\n", "\n", "def on_reset(b):\n", " global simplicits_state_0, simplicits_state_1, mpm_state_0, mpm_state_1\n", " sim_running[0] = False\n", " if sim_thread[0] is not None:\n", " sim_thread[0].join()\n", " simplicits_model.simplicits_scene.reset_scene()\n", " simplicits_state_0 = simplicits_model.state()\n", " simplicits_state_1 = simplicits_model.state()\n", " mpm_state_0 = mpm_model.state()\n", " mpm_state_1 = mpm_model.state()\n", " update_vis_mesh()\n", "\n", "\n", "buttons = [Button(description=x) for x in [\"Play\", \"Stop\", \"Reset\"]]\n", "buttons[0].on_click(on_play)\n", "buttons[1].on_click(on_stop)\n", "buttons[2].on_click(on_reset)\n", "\n", "update_vis_mesh()\n", "display(VBox(buttons), out)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.12.3" } }, "nbformat": 4, "nbformat_minor": 5 }