Skip to content
Skild AISkild AISteel

Nut bolt hydro

Mesh and SDF contacts

A nut on a bolt with hydroelastic mesh contacts.

0.00 s
The steel preset (the default), recordedon an NVIDIA RTX PRO 6000 Blackwell.

Run it

Terminal window
cd examples
# the default preset (steel), then the others
uv run shelf play nut_bolt_hydro
uv run shelf play nut_bolt_hydro --solver xpbd
uv run shelf play nut_bolt_hydro --solver mujoco
uv run shelf play nut_bolt_hydro --solver kamino
# change solver settings live
uv run shelf tune nut_bolt_hydro

Source: examples/shelf/03_mesh_sdf_contacts/nut_bolt_hydro/ · ported from newton/examples/contacts/example_nut_bolt_hydro.py @ skild-ai/newton 7565a26b1

Presets

PresetSolverStatusChecked
steeldefaultsteelworks2026-11-04, newton 7565a26b1, NVIDIA GeForce RTX 3090
xpbdxpbdworks2026-10-28, newton 7565a26b1, NVIDIA RTX A6000
mujocomujocoworks2026-10-28, newton 7565a26b1, NVIDIA RTX A6000
kaminokaminoworks2026-10-28, newton 7565a26b1, NVIDIA RTX A6000

Solvers with no row are untested on this example. The full settings are inpresets.toml.

Code

examples/shelf/03_mesh_sdf_contacts/nut_bolt_hydro/example.py
# SPDX-FileCopyrightText: Copyright (c) 2026 The Newton Developers
# SPDX-License-Identifier: Apache-2.0
# Modified for Steel Shelf: preset-driven execution and shared runner integration.


###########################################################################
# Example SDF Mesh Collision
#
# Demonstrates nut/bolt mesh collision using hydroelastic contacts.
#
# Command: python -m shelf run nut_bolt_hydro
#          python -m shelf run nut_bolt_hydro --solver kamino
#
###########################################################################


import argparse
import tempfile
from pathlib import Path


import numpy as np
import trimesh
import warp as wp


import newton
import newton.examples
from common import presets
from newton.geometry import HydroelasticSDF


# Assembly type for the nut and bolt
ASSEMBLY_STR = "m20_loose"


ISAACGYM_ENVS_REPO_URL = "https://github.com/isaac-sim/IsaacGymEnvs.git"
ISAACGYM_NUT_BOLT_FOLDER = "assets/factory/mesh/factory_nut_bolt"


SDF_MAX_RESOLUTION = 128
SDF_NARROW_BAND_RANGE = (-0.005, 0.005)
SDF_CONSTRUCTION_PADDING = 0.005
# Persist cooked SDFs across runs so the (slow) cook only happens once.
# Entries are content-addressed, so leftovers from older runs are harmless.
MESH_SDF_CACHE_DIR = Path(tempfile.gettempdir()) / "newton_sdf_cache"


SHAPE_CFG = newton.ModelBuilder.ShapeConfig(
    margin=0.0,
    mu=0.01,
    # Hydroelastic supplies the per-contact stiffness for the nut/bolt pair, so
    # ``ke``/``kd`` reach only the non-hydroelastic contacts. At the 1e10 ``kh``
    # default the thread contacts resolve over a ~95 ms solref time constant --
    # roughly 45 substeps -- and the nut sits visibly cocked under MuJoCo. XPBD
    # projects positions and is insensitive to this either way.
    kh=1e11,  # Hydroelastic contact stiffness
    ke=1e7,
    kd=1e4,
    gap=0.0,
    density=8000.0,
    mu_torsional=0.0,
    mu_rolling=0.0,
    is_hydroelastic=True,
)




# Demonstrate the user-facing pressure-callback API for the hydroelastic
# solver. The contact patch is defined as the iso-pressure surface
# ``p_a == p_b``; users supply a Warp ``@wp.func`` that maps a signed depth
# to a pressure value, plus a ``@wp.struct`` carrying any per-shape state it
# needs. The callback must be finite for any signed depth (positive or
# negative) and monotone non-increasing in ``signed_depth`` so the marching-
# cubes interpolation stays continuous across the patch boundary. Do not clip
# the non-contact side to zero pressure; with different shape stiffnesses, the
# iso-pressure surface can pass through that thin outside region.
#
# Here we re-implement the default linear law ``pressure = -kh * signed_depth``
# explicitly so the example exercises the user pathway. Nonlinear laws should
# similarly extend into ``signed_depth >= 0`` instead of flattening there.
@wp.struct
class LinearPressureData:
    shape_kh: wp.array[wp.float32]




@wp.func
def linear_pressure(signed_depth: wp.float32, shape_idx: wp.int32, data: LinearPressureData) -> wp.float32:
    return -data.shape_kh[shape_idx] * signed_depth




def add_mesh_object(
    builder: newton.ModelBuilder,
    mesh: newton.Mesh,
    transform: wp.transform,
    shape_cfg: newton.ModelBuilder.ShapeConfig | None = None,
    key: str | None = None,
    center_vec: wp.vec3 | None = None,
    scale: float = 1.0,
) -> int:
    """Add a mesh shape on a new body.


    Args:
        builder: Model builder that receives the body and shape.
        mesh: Mesh geometry with SDF data.
        transform: Body transform with position [m] and orientation.
        shape_cfg: Optional shape configuration.
        key: Optional body label.
        center_vec: Optional mesh center offset in local coordinates [m].
        scale: Uniform mesh scale [unitless].


    Returns:
        Created body index.
    """
    if center_vec is not None:
        center_world = wp.quat_rotate(transform.q, center_vec)
        transform = wp.transform(transform.p + center_world, transform.q)


    body = builder.add_body(label=key, xform=transform)
    builder.add_shape_mesh(body, mesh=mesh, scale=(scale, scale, scale), cfg=shape_cfg)
    return body




def load_mesh_with_sdf(
    mesh_file: str,
    scale: float = 1.0,
    center_origin: bool = True,
) -> tuple[newton.Mesh, wp.vec3]:
    """Load a triangle mesh and build an SDF.


    Args:
        mesh_file: Mesh file path.
        scale: Uniform mesh scale [unitless].
        center_origin: Whether to recenter mesh vertices about the AABB center.


    Returns:
        Tuple of ``(mesh, center_vec)`` where ``center_vec`` is the recenter offset [m].
    """
    mesh_data = trimesh.load(mesh_file, force="mesh")
    vertices = np.array(mesh_data.vertices, dtype=np.float32)
    indices = np.array(mesh_data.faces.flatten(), dtype=np.int32)
    center_vec = wp.vec3(0.0, 0.0, 0.0)


    if center_origin:
        min_extent = vertices.min(axis=0)
        max_extent = vertices.max(axis=0)
        center = (min_extent + max_extent) / 2
        vertices = vertices - center
        center_vec = wp.vec3(center) * float(scale)


    mesh = newton.Mesh(vertices, indices)
    mesh.build_sdf(
        max_resolution=SDF_MAX_RESOLUTION,
        narrow_band_range=SDF_NARROW_BAND_RANGE,
        margin=SDF_CONSTRUCTION_PADDING,
        scale=(scale, scale, scale),
        cache_dir=MESH_SDF_CACHE_DIR,
    )
    return mesh, center_vec




class Example:
    def __init__(self, viewer, args):
        self.preset = presets.load(__file__)
        newton.use_coord_layout_targets = True
        self.fps = self.preset.fps
        self.frame_dt = 1.0 / self.fps
        self.sim_time = 0.0
        self.sim_substeps = self.preset.substeps
        self.collide_every = self.preset.sim.need("collide_every")  # re-collide every K substeps
        self.sim_dt = self.frame_dt / self.sim_substeps


        self.world_count = args.world_count
        if self.world_count is None:
            self.world_count = 8 if self.preset.name == "kamino" else 20
        self.viewer = viewer
        self.solver_type = self.preset.name
        self.test_mode = args.test
        self.deterministic = args.deterministic
        self.deterministic_solver = args.deterministic_solver


        # Scene scaling factor (1.0 = original size)
        self.scene_scale = 1.0


        # Ground plane offset (negative = below origin)
        self.ground_plane_offset = -0.01


        # Grid dimensions for nut/bolt scene (number of assemblies in X and Y)
        self.num_per_world = args.num_per_world
        self.grid_x = int(np.ceil(np.sqrt(self.num_per_world)))
        self.grid_y = int(np.ceil(self.num_per_world / self.grid_x))


        # Contact budget per world. MuJoCo's njmax/nconmax are per-world limits
        # while CollisionPipeline's rigid_contact_max covers every world, so both
        # derive from this rather than from a fixed whole-scene pool. A threading
        # assembly peaks near 115 contacts, so this leaves ~9x headroom.
        self.contacts_per_world = 1024 * self.num_per_world


        # Maximum number of rigid contacts to allocate (limits memory usage)
        self.rigid_contact_max = self.contacts_per_world * self.world_count


        # Broad phase mode: NXN (O(N²)), SAP (O(N log N)), EXPLICIT (precomputed pairs)
        self.broad_phase_mode = "sap"


        world_builder = self._build_nut_bolt_scene()


        main_scene = newton.ModelBuilder()
        if self.solver_type == "kamino":
            newton.solvers.SolverKamino.register_custom_attributes(main_scene)
        main_scene.default_shape_cfg.gap = 0.001 * self.scene_scale
        # Add ground plane at z = ground_plane_offset.
        # For plane equation n·x + d = 0, with n=(0,0,1): z + d = 0, so z = -d.
        # Therefore d is the negative offset, and z = offset uses d = -offset.
        main_scene.add_shape_plane(
            plane=(0.0, 0.0, 1.0, -self.ground_plane_offset),
            width=0.0,
            length=0.0,
            label="ground_plane",
        )
        main_scene.replicate(world_builder, world_count=self.world_count)


        self.model = main_scene.finalize()


        # Configure the hydroelastic pipeline with our custom (still linear)
        # pressure callback. ``shape_kh`` reuses the per-shape stiffness already
        # stored on the model. Disable the marching-cubes edge clamp because
        # threading dynamics on the M20 helix are sensitive to the contact-
        # surface vertex bias the clamp introduces.
        pressure_data = LinearPressureData()
        pressure_data.shape_kh = self.model.shape_material_kh
        sdf_hydroelastic_config = HydroelasticSDF.Config(
            pressure_func=linear_pressure,
            pressure_data=pressure_data,
            mc_edge_clamp_min=0.0,
        )


        self.collision_pipeline = newton.CollisionPipeline(
            self.model,
            reduce_contacts=True,
            rigid_contact_max=self.rigid_contact_max,
            broad_phase=self.broad_phase_mode,
            sdf_hydroelastic_config=sdf_hydroelastic_config,
            deterministic=self.deterministic,
        )


        self.build_solver()


        self.state_0 = self.model.state()
        self.state_1 = self.model.state()
        self.control = self.model.control()


        newton.eval_fk(self.model, self.model.joint_q, self.model.joint_qd, self.state_0)


        self.contacts = self.collision_pipeline.contacts()
        self.collision_pipeline.collide(self.state_0, self.contacts)


        self.viewer.set_model(self.model)


        offset = 0.15 * self.scene_scale
        self.viewer.set_world_offsets((offset, offset, 0.0))
        camera_offset = np.sqrt(self.world_count) * offset * 1.25
        self.viewer.set_camera(pos=wp.vec3(camera_offset, -camera_offset, 0.5 * camera_offset), pitch=-15.0, yaw=135.0)


        # Initialize test tracking data (only in test mode for nut_bolt scene)
        self._init_test_tracking()


        self.capture()


    def build_solver(self):
        """Build the solver and everything bound to its instance; the tuner calls this again to swap it."""
        # Create solver based on user choice
        if self.solver_type == "xpbd":
            self.solver = self.preset.make_solver(
                self.model,
                deterministic=wp.DeterministicMode.RUN_TO_RUN
                if self.deterministic_solver
                else wp.DeterministicMode.NOT_GUARANTEED,
            )
        elif self.solver_type == "mujoco":
            self.solver = self.preset.make_solver(
                self.model,
                njmax=self.contacts_per_world,
                nconmax=self.contacts_per_world,
                deterministic=wp.DeterministicMode.RUN_TO_RUN
                if self.deterministic_solver
                else wp.DeterministicMode.NOT_GUARANTEED,
            )
        elif self.solver_type == "kamino":
            self.solver = self.preset.make_solver(self.model)
        elif self.solver_type == "steel":
            # Free-body contact rows sized per assembly; the contact law comes from the preset.
            self.solver = self.preset.make_solver(
                self.model,
                free_body_constraint_max=3 * 256 * self.num_per_world,
            )
        else:
            raise ValueError(f"Unknown solver '{self.solver_type}'")


    def _build_nut_bolt_scene(self) -> newton.ModelBuilder:
        print("Downloading nut/bolt assets...")
        asset_path = newton.examples.download_external_git_folder(ISAACGYM_ENVS_REPO_URL, ISAACGYM_NUT_BOLT_FOLDER, ref="aeed298638a1f7b5421b38f5f3cc2d1079b6d9c3")
        print(f"Assets downloaded to: {asset_path}")


        world_builder = newton.ModelBuilder()
        world_builder.default_shape_cfg.gap = 0.001 * self.scene_scale


        bolt_file = str(asset_path / f"factory_bolt_{ASSEMBLY_STR}.obj")
        nut_file = str(asset_path / f"factory_nut_{ASSEMBLY_STR}_subdiv_3x.obj")
        bolt_mesh, bolt_center = load_mesh_with_sdf(bolt_file, scale=self.scene_scale, center_origin=True)
        nut_mesh, nut_center = load_mesh_with_sdf(nut_file, scale=self.scene_scale, center_origin=True)


        # Spacing between assemblies in the grid
        spacing = 0.1 * self.scene_scale


        # Create grid of nut/bolt assemblies
        count = 0
        for i in range(self.grid_x):
            if count >= self.num_per_world:
                break
            for j in range(self.grid_y):
                if count >= self.num_per_world:
                    break
                # Center the grid around origin
                x_offset = (i - (self.grid_x - 1) / 2.0) * spacing
                y_offset = (j - (self.grid_y - 1) / 2.0) * spacing


                # Add bolt at grid position
                bolt_xform = wp.transform(wp.vec3(x_offset, y_offset, 0.0 * self.scene_scale), wp.quat_identity())
                add_mesh_object(
                    world_builder,
                    bolt_mesh,
                    bolt_xform,
                    SHAPE_CFG,
                    key=f"bolt_{i}_{j}",
                    center_vec=bolt_center,
                    scale=self.scene_scale,
                )


                # Add nut above bolt at grid position
                nut_xform = wp.transform(
                    wp.vec3(x_offset, y_offset, 0.041 * self.scene_scale),
                    wp.quat_from_axis_angle(wp.vec3(0.0, 0.0, 1.0), np.pi / 8),
                )
                add_mesh_object(
                    world_builder,
                    nut_mesh,
                    nut_xform,
                    SHAPE_CFG,
                    key=f"nut_{i}_{j}",
                    center_vec=nut_center,
                    scale=self.scene_scale,
                )
                count += 1


        return world_builder


    def capture(self):
        with wp.ScopedCapture() as capture:
            self.simulate()
        self.graph = capture.graph


    def simulate(self):
        for sub in range(self.sim_substeps):
            # Refresh contacts every K substeps so contact normals stay
            # aligned with the threading rotation.
            if sub % self.collide_every == 0:
                self.collision_pipeline.collide(self.state_0, self.contacts)
            self.state_0.clear_forces()


            self.viewer.apply_forces(self.state_0)
            self.solver.step(self.state_0, self.state_1, self.control, self.contacts, self.sim_dt)


            if self.sim_substeps % 2 == 1 and sub == self.sim_substeps - 1:
                # an odd count under graph capture: copy the last result into state_0 instead of swapping,
                # so every replay of the captured graph starts from the newest state
                self.state_0.assign(self.state_1)
            else:
                self.state_0, self.state_1 = self.state_1, self.state_0


    def step(self):
        if self.graph:
            wp.capture_launch(self.graph)
        else:
            self.simulate()
        if self.solver_type == "steel":
            self.solver.validate_contact_compliance()


        self.sim_time += self.frame_dt


        # Track transforms for test validation
        self._track_test_data()


    def render(self):
        self.viewer.begin_frame(self.sim_time)
        self.viewer.log_state(self.state_0)
        self.viewer.log_contacts(self.contacts, self.state_0)
        self.viewer.end_frame()


    def _init_test_tracking(self):
        """Initialize tracking data for test validation."""
        if not self.test_mode:
            self.bolt_body_indices = None
            self.nut_body_indices = None
            return


        # Find bolt and nut body indices by key
        self.bolt_body_indices = []
        self.nut_body_indices = []


        for i in range(self.grid_x):
            for j in range(self.grid_y):
                bolt_key = f"bolt_{i}_{j}"
                nut_key = f"nut_{i}_{j}"


                if bolt_key in self.model.body_label:
                    self.bolt_body_indices.append(self.model.body_label.index(bolt_key))
                if nut_key in self.model.body_label:
                    self.nut_body_indices.append(self.model.body_label.index(nut_key))


        # Store initial transforms
        body_q = self.state_0.body_q.numpy()
        self.bolt_initial_transforms = [body_q[idx].copy() for idx in self.bolt_body_indices]
        self.nut_initial_transforms = [body_q[idx].copy() for idx in self.nut_body_indices]


        # Track maximum rotation change and z displacement for nuts
        self.nut_max_rotation_change = [0.0] * len(self.nut_body_indices)
        self.nut_min_z = [body_q[idx][2] for idx in self.nut_body_indices]


    def _track_test_data(self):
        """Track transforms for test validation (called each step in test mode)."""
        if not self.test_mode:
            return


        body_q = self.state_0.body_q.numpy()


        # Track nut rotation and z position
        for i, nut_idx in enumerate(self.nut_body_indices):
            current_q = body_q[nut_idx]
            initial_q = self.nut_initial_transforms[i]


            # Compute rotation change using quaternion dot product
            # |q1 · q2| = cos(theta/2), where theta is the angle between orientations
            q_current = current_q[3:7]  # quaternion part (x, y, z, w)
            q_initial = initial_q[3:7]
            dot = abs(np.dot(q_current, q_initial))
            dot = min(dot, 1.0)  # Clamp for numerical stability
            rotation_angle = 2.0 * np.arccos(dot)
            self.nut_max_rotation_change[i] = max(self.nut_max_rotation_change[i], rotation_angle)


            # Track minimum z (nuts should move down)
            self.nut_min_z[i] = min(self.nut_min_z[i], current_q[2])


    def test_final(self):
        """Verify simulation state after example completes.


        - Bolts should stay approximately in place (limited displacement)
        - Nuts should rotate (thread engagement) and move slightly downward
        """
        body_q = self.state_0.body_q.numpy()


        # Check bolts stayed in place
        max_bolt_displacement = 0.02  # 2 cm tolerance
        for i, bolt_idx in enumerate(self.bolt_body_indices):
            current_pos = body_q[bolt_idx][:3]
            initial_pos = self.bolt_initial_transforms[i][:3]
            displacement = np.linalg.norm(current_pos - initial_pos)
            assert displacement < max_bolt_displacement, (
                f"Bolt {i}: displaced too much. "
                f"Displacement={displacement:.4f} (max allowed={max_bolt_displacement:.4f})"
            )


        # The 45 degree threshold catches stalled thread engagement while
        # allowing observed solver/contact-count variation.
        min_rotation_threshold = np.radians(45.0)
        min_descent = 0.005
        for i in range(len(self.nut_body_indices)):
            # Check rotation occurred
            max_rotation = self.nut_max_rotation_change[i]
            assert max_rotation > min_rotation_threshold, (
                f"Nut {i}: did not rotate enough. "
                f"Max rotation={np.degrees(max_rotation):.2f} degrees "
                f"(expected > {np.degrees(min_rotation_threshold):.2f} degrees)"
            )


            # Check nut moved downward (min_z should be less than initial z)
            initial_z = self.nut_initial_transforms[i][2]
            min_z = self.nut_min_z[i]
            descent = initial_z - min_z
            assert descent > min_descent, (
                f"Nut {i}: did not move down enough. Descent={descent:.4f} (expected > {min_descent:.4f})"
            )


    @staticmethod
    def create_parser():
        parser = newton.examples.create_parser()
        newton.examples.add_world_count_arg(parser)
        parser.set_defaults(world_count=None)
        parser.add_argument(
            "--deterministic",
            action=argparse.BooleanOptionalAction,
            default=False,
            help=(
                "Make contact generation and ordering reproducible across runs on the same GPU. "
                "Costs a few percent of step time. Needs a world count small enough to keep the "
                "hydroelastic face-contact buffer under the 2^20 deterministic contact-id limit."
            ),
        )
        parser.add_argument(
            "--deterministic-solver",
            action=argparse.BooleanOptionalAction,
            default=False,
            help=(
                "Additionally make the solver bit-exact. Separate from --deterministic because "
                "it instruments every atomic scatter in the solver and costs ~7x step time, "
                "while contact ordering is what varies between runs."
            ),
        )
        parser.add_argument(
            "--num-per-world",
            type=int,
            default=1,
            help="Number of assemblies per world.",
        )
        return parser




if __name__ == "__main__":
    parser = Example.create_parser()
    viewer, args = newton.examples.init(parser)


    newton.examples.run(Example(viewer, args), args)

example.py on GitHub

Steel docs and tutorials are licensed under Apache-2.0. Third-party notices. "Skild", "Skild AI" and the Skild logo are trademarks of Skild AI and are not licensed under Apache-2.0.