Basic dzhanibekov

Basics

The Dzhanibekov effect on a free T-shaped rigid body.

0.00 s

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

Run it

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

Source: examples/shelf/00_basics/basic_dzhanibekov/ · ported from newton/examples/basic/example_basic_dzhanibekov.py @ skild-ai/newton 7565a26b1

Presets

Preset	Solver	Status	Checked
`steel`default	steel	works	2026-11-04, newton 7565a26b1, NVIDIA GeForce RTX 3090
`xpbd`	xpbd	works	2026-10-28, newton 7565a26b1, NVIDIA RTX A6000
`vbd`	vbd	works	2026-10-28, newton 7565a26b1, NVIDIA RTX A6000
`mujoco`	mujoco	works	2026-10-28, newton 7565a26b1, NVIDIA RTX A6000
`kamino`	kamino	works	2026-10-28, newton 7565a26b1, NVIDIA RTX A6000

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

# 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 Basic Dzhanibekov
#
# Demonstrates the Dzhanibekov effect on a free T-shaped rigid body.
# The body is composed of a horizontal bar (一) and a vertical stem (|).
# The initial angular velocity is about the stem axis, gravity is
# disabled, and the body is advanced by a rigid-body solver.
#
# Command: python -m shelf run basic_dzhanibekov
# XPBD: python -m shelf run basic_dzhanibekov --solver xpbd
# MuJoCo: python -m shelf run basic_dzhanibekov --solver mujoco
# Kamino: python -m shelf run basic_dzhanibekov --solver kamino
#
###########################################################################

import numpy as np
import warp as wp

import newton
import newton.examples
from common import presets

_STEM_RADIUS = 0.10
_STEM_HALF_HEIGHT = 0.5
_STEM_DENSITY = 50.0

_BAR_RADIUS = 0.20
_BAR_HALF_HEIGHT = 1.0
_BAR_DENSITY = 100.0


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.sim_dt = self.frame_dt / self.sim_substeps

        self.viewer = viewer
        self.solver_type = self.preset.name
        self.min_stem_axis_y = 1.0

        builder = newton.ModelBuilder(up_axis=newton.Axis.Z)
        if self.solver_type == "kamino":
            newton.solvers.SolverKamino.register_custom_attributes(builder)

        self.body = builder.add_body(
            xform=wp.transform(wp.vec3(0.0, 0.0, 2.0), wp.quat_identity()),
            label="dzhanibekov_t",
        )
        self.free_joint = len(builder.joint_type) - 1

        # Build T-shaped rigid body: bar (一) + stem (|)
        # Cylinders default to Z axis; rotate bar to X and stem to Y.
        q_bar = wp.quat_from_axis_angle(wp.vec3(0.0, 1.0, 0.0), -wp.pi / 2.0)  # Z → X
        q_stem = wp.quat_from_axis_angle(wp.vec3(1.0, 0.0, 0.0), wp.pi / 2.0)  # Z → Y

        bar_center = np.array((0.0, 0.0, 0.0), dtype=np.float32)
        stem_center = np.array((0.0, -(_BAR_RADIUS + _STEM_HALF_HEIGHT), 0.0), dtype=np.float32)
        bar_mass = 2.0 * wp.pi * _BAR_RADIUS**2 * _BAR_HALF_HEIGHT * _BAR_DENSITY
        stem_mass = 2.0 * wp.pi * _STEM_RADIUS**2 * _STEM_HALF_HEIGHT * _STEM_DENSITY
        center_of_mass = (bar_center * bar_mass + stem_center * stem_mass) / (bar_mass + stem_mass)

        cfg_bar = builder.default_shape_cfg.copy()
        cfg_bar.density = _BAR_DENSITY
        cfg_bar.has_shape_collision = False
        builder.add_shape_cylinder(
            self.body,
            xform=wp.transform(wp.vec3(*(bar_center - center_of_mass)), q_bar),
            radius=_BAR_RADIUS,
            half_height=_BAR_HALF_HEIGHT,
            cfg=cfg_bar,
            color=wp.vec3(0.75, 0.18, 0.12),
        )

        cfg_stem = builder.default_shape_cfg.copy()
        cfg_stem.density = _STEM_DENSITY
        cfg_stem.has_shape_collision = False
        builder.add_shape_cylinder(
            self.body,
            xform=wp.transform(wp.vec3(*(stem_center - center_of_mass)), q_stem),
            radius=_STEM_RADIUS,
            half_height=_STEM_HALF_HEIGHT,
            cfg=cfg_stem,
            color=wp.vec3(0.15, 0.42, 0.78),
        )

        builder.body_qd[self.body] = wp.spatial_vector(0.0, 0.0, 0.0, 0.0, 20.0, 0.0)
        builder.joint_qd = np.array(builder.body_qd).flatten().tolist()

        if self.solver_type == "vbd":
            builder.color()

        builder.add_ground_plane()

        self.model = builder.finalize()
        self.model.set_gravity((0.0, 0.0, 0.0))

        self.build_solver()

        self.state_0 = self.model.state()
        self.state_1 = self.model.state()
        self.control = self.model.control()
        self.initial_rotational_energy = self._rotational_energy() if self.solver_type == "kamino" else None

        self.initial_body_q = self.state_0.body_q.numpy().copy()
        self.body_com = self.model.body_com.numpy()[self.body].copy()
        self.initial_com = (
            self.initial_body_q[self.body][:3]
            + np.array(wp.quat_to_matrix(wp.quat(*self.initial_body_q[self.body][3:7])), dtype=np.float32).reshape(3, 3)
            @ self.body_com
        )

        self.viewer.set_model(self.model)

        self.capture()

    def build_solver(self):
        """Build the solver and everything bound to its instance; the tuner calls this again to swap it."""
        if self.solver_type == "vbd":
            self.solver = self.preset.make_solver(self.model)
            self.collision_pipeline = newton.CollisionPipeline(self.model)
            self.contacts = self.collision_pipeline.contacts()
        elif self.solver_type in ("xpbd", "steel"):
            self.solver = self.preset.make_solver(self.model)
            self.collision_pipeline = newton.CollisionPipeline(self.model)
            self.contacts = self.collision_pipeline.contacts()
        elif self.solver_type == "mujoco":
            self.solver = self.preset.make_solver(self.model)
            self.collision_pipeline = None
            self.contacts = None
        elif self.solver_type == "kamino":
            self.solver = self.preset.make_solver(self.model)
            self.collision_pipeline = None
            self.contacts = None
        else:
            raise ValueError(
                f"Unknown solver type: {self.solver_type}. Choose from 'vbd', 'xpbd', 'steel', 'mujoco', or 'kamino'."
            )

    def capture(self):
        if wp.get_device().is_cuda:
            with wp.ScopedCapture() as capture:
                self.simulate()
            self.graph = capture.graph
        else:
            self.graph = None

    def simulate(self):
        for substep in range(self.sim_substeps):
            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 substep == 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()

        self.sim_time += self.frame_dt
        body_q = self.state_0.body_q.numpy()
        local_y_world = np.array(wp.quat_to_matrix(wp.quat(*body_q[self.body][3:7])), dtype=np.float32).reshape(3, 3)[
            :, 1
        ]
        self.min_stem_axis_y = min(self.min_stem_axis_y, float(local_y_world[1]))

    def _rotational_energy(self):
        body_q = self.state_0.body_q.numpy()[self.body]
        omega = self.state_0.body_qd.numpy()[self.body, 3:6]
        rotation = np.array(wp.quat_to_matrix(wp.quat(*body_q[3:7])), dtype=np.float32).reshape(3, 3)
        inertia = self.model.body_inertia.numpy()[self.body]
        return float(0.5 * omega @ (rotation @ inertia @ rotation.T) @ omega)

    def test_final(self):
        """Verify the free body flips without gaining rotational energy."""
        if self.initial_rotational_energy is not None:
            energy_ratio = self._rotational_energy() / self.initial_rotational_energy
            if not np.isfinite(energy_ratio) or abs(energy_ratio - 1.0) > 0.01:
                raise ValueError(f"Torque-free body changed rotational energy by {energy_ratio - 1.0:.1%}")
        body_q = self.state_0.body_q.numpy()
        body_qd = self.state_0.body_qd.numpy()
        if not np.isfinite(body_q[self.body]).all() or not np.isfinite(body_qd[self.body]).all():
            raise ValueError("Dzhanibekov rigid body state contains non-finite values")
        final_com = (
            body_q[self.body][:3]
            + np.array(wp.quat_to_matrix(wp.quat(*body_q[self.body][3:7])), dtype=np.float32).reshape(3, 3)
            @ self.body_com
        )
        if np.linalg.norm(final_com - self.initial_com) >= 1.0e-3:
            raise ValueError("Free rigid body developed linear drift")
        if self.min_stem_axis_y >= 0.0:
            raise ValueError("Dzhanibekov flip did not occur")

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


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

example.py on GitHub

example = "basic_dzhanibekov"
source = "newton/examples/basic/example_basic_dzhanibekov.py @ skild-ai/newton 7565a26b1"
default = "steel"

[check]
# Arguments `python -m shelf check` passes in test mode (from Newton's test_examples.py).
args = []

[gallery]
short = "The Dzhanibekov effect on a free T-shaped rigid body."

[xpbd]
status = "works"
sim = { fps = 60, substeps = 4 }
solver = { iterations = 10 }
checked = "2026-10-28, newton 7565a26b1, NVIDIA RTX A6000"

[vbd]
status = "works"
sim = { fps = 60, substeps = 4 }
solver = { iterations = 4 }
checked = "2026-10-28, newton 7565a26b1, NVIDIA RTX A6000"

[mujoco]
status = "works"
sim = { fps = 60, substeps = 4 }
solver = { disable_contacts = true, solver = "newton", integrator = "implicitfast", iterations = 10 }
checked = "2026-10-28, newton 7565a26b1, NVIDIA RTX A6000"

[featherstone]
status = "untested"

[semi_implicit]
status = "untested"

[feather_pgs]
status = "untested"

[steel]
status = "works"
# 4 substeps, immediate response: 2 substeps (0.30 ms) delays the flip past the 100 frames the test allows, and
# odd counts freeze this example's graph; the body has no contact rows, so iterations cost nothing, and
# propagation_colored, split mode or serial streams are all slower here.
sim = { fps = 60, substeps = 4 }
solver = { solve_mode = "matrix_free", iterations = 4, articulated_contact_response = "immediate" }
checked = "2026-11-04, newton 7565a26b1, NVIDIA GeForce RTX 3090"

[kamino]
status = "works"
sim = { fps = 60, substeps = 4 }
solver = { dynamics_solver = "dvi", sparse_dynamics = true, sparse_jacobian = true }
checked = "2026-10-28, newton 7565a26b1, NVIDIA RTX A6000"

[implicit_mpm]
status = "untested"

[style3d]
status = "untested"

presets.toml 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.