diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md
index 8feaf9cd..6ede5d6c 100644
--- a/CONTRIBUTING.md
+++ b/CONTRIBUTING.md
@@ -7,9 +7,10 @@ Thank you for your interest in contributing to Genesis! We welcome contributions
 - [Submitting Code Changes](#submitting-code-changes)
 - [Reviewing and Merging](#reviewing-and-merging)
 - [Questions and Discussions](#questions-and-discussions)
+
 ---
 
-### Reporting Bugs
+## Reporting Bugs
 
 - Before reporting a bug, please search through existing issues to check if it has already been reported.
 
@@ -36,36 +37,44 @@ Thank you for your interest in contributing to Genesis! We welcome contributions
   Add any other context about the problem here.
   ```
 
-### Suggesting Features
+## Suggesting Features
+
 - If you have a feature idea, please create an issue labeled `enhancement`.
 - In the created issue, please provide context, expected outcomes, and potential.
 
-### Submitting Code Changes
+## Submitting Code Changes
 
 - We use the `pre-commit` configuration to automatically clean up code before committing. Install and run `pre-commit` as follows:
   1. Install `pre-commit`:
+
      ```bash
      pip install pre-commit
      ```
+
   2. Install hooks from the configuration file:
+
      ```bash
      pre-commit install
      ```
+
      After this, `pre-commit` will automatically check and clean up code whenever you make a commit.
 - (Optional) You can run CI tests locally to ensure you pass the online CI checks.
+
   ```python
   python -m unittest discover tests
   ```
+
 - In the title of your Pull Request, please include [BUG FIX], [FEATURE] or [MISC] to indicate the purpose.
 - In the description, please provide example code or commands for testing.
 
-### Reviewing and Merging
+## Reviewing and Merging
 
 - PRs require at least one approval before merging.
 - Automated checks (e.g., CI tests) must pass.
 - Use `Squash and Merge` for a clean commit history.
 
-### Questions and Discussions
+## Questions and Discussions
+
 - Use [Discussions](https://github.com/Genesis-Embodied-AI/Genesis/discussions) for open-ended topics.
 <!-- 
 ### Join Us
diff --git a/README.md b/README.md
index 25c5c5a2..47ac6e54 100644
--- a/README.md
+++ b/README.md
@@ -9,43 +9,56 @@
 
 [![README in English](https://img.shields.io/badge/English-d9d9d9)](./README.md)
 [![简体中文版自述文件](https://img.shields.io/badge/简体中文-d9d9d9)](./README_CN.md)
+[![日本語版 README](https://img.shields.io/badge/日本語-d9d9d9)](./README_JA.md)
 
-# What is Genesis?
+# Genesis
+## 🔥 News
+- [2024-12-25] Added a [docker](#docker) including support for the ray-tracing renderer
+- [2024-12-24] Added guidelines for [contributing to Genesis](https://github.com/Genesis-Embodied-AI/Genesis/blob/main/CONTRIBUTING.md)
 
-Genesis is a physics platform designed for general purpose *Robotics/Embodied AI/Physical AI* applications. It is simultaneously multiple things:
+## Table of Contents
+
+1. [What is Genesis?](#what-is-genesis)
+2. [Key Features](#key-features)
+3. [Quick Installation](#quick-installation)
+4. [Docker](#docker)
+5. [Documentation](#documentation)
+6. [Contributing to Genesis](#contributing-to-genesis)
+7. [Support](#support)
+8. [License and Acknowledgments](#license-and-acknowledgments)
+9. [Associated Papers](#associated-papers)
+10. [Citation](#citation)
+
+## What is Genesis?
+
+Genesis is a physics platform designed for general-purpose *Robotics/Embodied AI/Physical AI* applications. It is simultaneously multiple things:
 
 1. A **universal physics engine** re-built from the ground up, capable of simulating a wide range of materials and physical phenomena.
 2. A **lightweight**, **ultra-fast**, **pythonic**, and **user-friendly** robotics simulation platform.
 3. A powerful and fast **photo-realistic rendering system**.
 4. A **generative data engine** that transforms user-prompted natural language description into various modalities of data.
 
-Powered by a universal physics engine re-designed and re-built from the ground up, Genesis integrates various physics solvers and their coupling into a unified framework. This core physics engine is further enhanced by a generative agent framework that operates at an upper level, aiming towards fully **automated data generation** for robotics and beyond.
+Genesis aims to:
 
-Currently, we are open-sourcing the **underlying physics engine and the simulation platform**. Our generative framework is a modular system that incorporates many different generative modules, each handling a certain range of data modalities, routed by a high level agent. Some of the modules integrated existing papers and some are still under submission. Access to our generative feature will be gradually rolled out in the near future. If you are interested, feel free to explore more the [paper list](#papers-behind-genesis) below.
-
-Genesis is built and will continuously evolve with the following ***long-term missions***:
-
-1. **Lowering the barrier** to using physics simulations and making robotics research accessible to everyone. (See our [commitment](https://genesis-world.readthedocs.io/en/latest/user_guide/overview/mission.html))
-2. **Unifying a wide spectrum of state-of-the-art physics solvers** into a single framework, allowing re-creating the whole physical world in a virtual realm with the highest possible physical, visual and sensory fidelity, using the most advanced simulation techniques.
-3. **Minimizing human effort** in collecting and generating data for robotics and other domains, letting the data flywheel spin on its own.
+- **Lower the barrier** to using physics simulations, making robotics research accessible to everyone. See our [mission statement](https://genesis-world.readthedocs.io/en/latest/user_guide/overview/mission.html).
+- **Unify diverse physics solvers** into a single framework to recreate the physical world with the highest fidelity.
+- **Automate data generation**, reducing human effort and letting the data flywheel spin on its own.
 
 Project Page: <https://genesis-embodied-ai.github.io/>
 
 ## Key Features
 
-- **Speed**: Genesis delivers an unprecedented simulation speed -- over 43 million FPS when simulating a Franka robotic arm with a single RTX 4090 (430,000 times faster than real-time).
-- **Cross-platform**: Genesis runs natively across different systems (Linux, MacOS, Windows), and across different compute backend (CPU, Nvidia GPU, AMD GPU, Apple Metal).
-- **Unification of various physics solvers**: Genesis develops a unified simulation framework that integrates various physics solvers: Rigid body, MPM, SPH, FEM, PBD, Stable Fluid.
-- **Support a wide range of material models**:  Genesis supports simulation (and the coupling) of rigid and articulated bodies, various types of liquids, gaseous phenomenon, deformable objects, thin-shell objects and granular materials.
-- **Support for a wide range of robots**: Robot arm, legged robot, drone, *soft robot*, etc., and extensive support for loading different file types: `MJCF (.xml)`, `URDF`, `.obj`, `.glb`, `.ply`, `.stl`, etc.
-- **Photorealistic and high-performance ray-tracer**: Genesis supports native ray-tracing based rendering.
-- **Differentiability**: Genesis is designed to be fully compatible with differentiable simulation. Currently, our MPM solver and Tool Solver are differentiable, and differentiability for other solvers will be added soon (starting with rigid-body simulation).
-- **Physics-based Tactile Sensor**: Genesis involves a physics-based and differentiable [tactile sensor simulation module](https://github.com/Genesis-Embodied-AI/DiffTactile). This will be integrated to the public version soon (expected in version 0.3.0).
-- **User-friendliness**: Genesis is designed in a way to make using simulation as simple as possible. From installation to API design, if there's anything you found counter-intuitive or difficult to use, please [let us know](https://github.com/Genesis-Embodied-AI/Genesis/issues).
-
-## Getting Started
+- **Speed**: Over 43 million FPS when simulating a Franka robotic arm with a single RTX 4090 (430,000 times faster than real-time).
+- **Cross-platform**: Runs on Linux, macOS, Windows, and supports multiple compute backends (CPU, Nvidia/AMD GPUs, Apple Metal).
+- **Integration of diverse physics solvers**: Rigid body, MPM, SPH, FEM, PBD, Stable Fluid.
+- **Wide range of material models**: Simulation and coupling of rigid bodies, liquids, gases, deformable objects, thin-shell objects, and granular materials.
+- **Compatibility with various robots**: Robotic arms, legged robots, drones, *soft robots*, and support for loading `MJCF (.xml)`, `URDF`, `.obj`, `.glb`, `.ply`, `.stl`, and more.
+- **Photo-realistic rendering**: Native ray-tracing-based rendering.
+- **Differentiability**: Genesis is designed to be fully differentiable. Currently, our MPM solver and Tool Solver support differentiability, with other solvers planned for future versions (starting with rigid & articulated body solver).
+- **Physics-based tactile simulation**: Differentiable [tactile sensor simulation](https://github.com/Genesis-Embodied-AI/DiffTactile) coming soon (expected in version 0.3.0).
+- **User-friendliness**: Designed for simplicity, with intuitive installation and APIs.
 
-### Quick Installation
+## Quick Installation
 
 Genesis is available via PyPI:
 
@@ -55,40 +68,69 @@ pip install genesis-world  # Requires Python >=3.9;
 
 You also need to install **PyTorch** following the [official instructions](https://pytorch.org/get-started/locally/).
 
-If you would like to try out the latest version, we suggest you to git clone from the repo and do `pip install -e .` instead of via PyPI.
+For the latest version, clone the repository and install locally:
 
-### Documentation
+```bash
+git clone https://github.com/Genesis-Embodied-AI/Genesis.git
+cd Genesis
+pip install -e .
+```
+
+## Docker
 
-Please refer to our [documentation site (English)](https://genesis-world.readthedocs.io/en/latest/user_guide/index.html) / [(Chinese)](https://genesis-world.readthedocs.io/zh-cn/latest/user_guide/index.html) for detailed installation steps, tutorials and API references. 
+If you want to use Genesis from Docker, you can first build the Docker image as:
+
+```bash
+docker build -t genesis -f docker/Dockerfile docker
+```
+
+Then you can run the examples inside the docker image (mounted to `/workspace/examples`):
+
+```bash
+xhost +local:root # Allow the container to access the display
+
+docker run --gpus all --rm -it \
+-e DISPLAY=$DISPLAY \
+-v /tmp/.X11-unix/:/tmp/.X11-unix \
+-v $PWD:/workspace \
+genesis
+```
+
+## Documentation
+
+Comprehensive documentation is available in [English](https://genesis-world.readthedocs.io/en/latest/user_guide/index.html) and [Chinese](https://genesis-world.readthedocs.io/zh-cn/latest/user_guide/index.html). This includes detailed installation steps, tutorials, and API references.
 
 ## Contributing to Genesis
 
-The goal of the Genesis project is to build a fully transparent, user-friendly ecosystem where contributors from both robotics and computer graphics can **come together to collaboratively create a high-efficiency, realistic (both physically and visually) virtual world for robotics research and beyond**.
+The Genesis project is an open and collaborative effort. We welcome all forms of contributions from the community, including:
 
-We sincerely welcome *any forms of contributions* from the community to make the world a better place for robots. From **pull requests** for new features, **bug reports**, to even tiny **suggestions** that will make Genesis API more intuitive, all are wholeheartedly appreciated!
+- **Pull requests** for new features or bug fixes.
+- **Bug reports** through GitHub Issues.
+- **Suggestions** to improve Genesis's usability.
 
-## Support
+Refer to our [contribution guide](https://github.com/Genesis-Embodied-AI/Genesis/blob/main/CONTRIBUTING.md) for more details.
 
-- Please use Github [Issues](https://github.com/Genesis-Embodied-AI/Genesis/issues) for bug reports and feature requests.
+## Support
 
-- Please use GitHub [Discussions](https://github.com/Genesis-Embodied-AI/Genesis/discussions) for discussing ideas, and asking questions.
+- Report bugs or request features via GitHub [Issues](https://github.com/Genesis-Embodied-AI/Genesis/issues).
+- Join discussions or ask questions on GitHub [Discussions](https://github.com/Genesis-Embodied-AI/Genesis/discussions).
 
-## License and Acknowledgment
+## License and Acknowledgments
 
 The Genesis source code is licensed under Apache 2.0.
-The development of Genesis won't be possible without these amazing open-source projects:
 
-- [Taichi](https://github.com/taichi-dev/taichi): for providing a high-performance cross-platform compute backend. Kudos to all the members providing technical support from taichi!
-- [FluidLab](https://github.com/zhouxian/FluidLab) for providing a reference MPM solver implementation
-- [SPH_Taichi](https://github.com/erizmr/SPH_Taichi) for providing a reference SPH solver implementation
-- [Ten Minute Physics](https://matthias-research.github.io/pages/tenMinutePhysics/index.html) and [PBF3D](https://github.com/WASD4959/PBF3D) for providing a reference PBD solver implementation
-- [MuJoCo](https://github.com/google-deepmind/mujoco) and [Brax](https://github.com/google/brax) for providing reference for rigid body dynamics
-- [libccd](https://github.com/danfis/libccd) for providing reference for collision detection
-- [PyRender](https://github.com/mmatl/pyrender) for rasterization-based renderer
-- [LuisaCompute](https://github.com/LuisaGroup/LuisaCompute) and [LuisaRender](https://github.com/LuisaGroup/LuisaRender) for its ray-tracing DSL
-- [trimesh](https://github.com/mikedh/trimesh), [PyMeshLab](https://github.com/cnr-isti-vclab/PyMeshLab) and [CoACD](https://github.com/SarahWeiii/CoACD) for geometry processing
+Genesis's development has been made possible thanks to these open-source projects:
 
-## Papers behind Genesis
+- [Taichi](https://github.com/taichi-dev/taichi): High-performance cross-platform compute backend. Kudos to the Taichi team for their technical support!
+- [FluidLab](https://github.com/zhouxian/FluidLab): Reference MPM solver implementation.
+- [SPH_Taichi](https://github.com/erizmr/SPH_Taichi): Reference SPH solver implementation.
+- [Ten Minute Physics](https://matthias-research.github.io/pages/tenMinutePhysics/index.html) and [PBF3D](https://github.com/WASD4959/PBF3D): Reference PBD solver implementations.
+- [MuJoCo](https://github.com/google-deepmind/mujoco): Reference for rigid body dynamics.
+- [libccd](https://github.com/danfis/libccd): Reference for collision detection.
+- [PyRender](https://github.com/mmatl/pyrender): Rasterization-based renderer.
+- [LuisaCompute](https://github.com/LuisaGroup/LuisaCompute) and [LuisaRender](https://github.com/LuisaGroup/LuisaRender): Ray-tracing DSL.
+
+## Associated Papers
 
 Genesis is a large scale effort that integrates state-of-the-art technologies of various existing and on-going research work into a single system. Here we include a non-exhaustive list of all the papers that contributed to the Genesis project in one way or another:
 
@@ -116,7 +158,7 @@ Genesis is a large scale effort that integrates state-of-the-art technologies of
 
 ## Citation
 
-If you used Genesis in your research, we would appreciate it if you could cite it. We are still working on a technical report, and before it's public, you could consider citing:
+If you use Genesis in your research, please consider citing:
 
 ```bibtex
 @software{Genesis,
@@ -126,4 +168,3 @@ If you used Genesis in your research, we would appreciate it if you could cite i
   year = {2024},
   url = {https://github.com/Genesis-Embodied-AI/Genesis}
 }
-```
diff --git a/README_CN.md b/README_CN.md
index 41e8c424..bb852175 100644
--- a/README_CN.md
+++ b/README_CN.md
@@ -1,3 +1,5 @@
+<div align="center">
+
 ![Genesis](imgs/big_text.png)
 
 ![Teaser](imgs/teaser.png)
@@ -9,27 +11,39 @@
 
 [![README in English](https://img.shields.io/badge/English-d9d9d9)](./README.md)
 [![简体中文版自述文件](https://img.shields.io/badge/简体中文-d9d9d9)](./README_CN.md)
+[![日本語版 README](https://img.shields.io/badge/日本語-d9d9d9)](./README_JA.md)
+
+</div>
 
-# 概述
+# Genesis 通用物理引擎
 
-Genesis 是一个为通用 *机器人/嵌入式 AI/物理 AI* 应用设计的物理平台。它同时具备多种功能：
+## 目录
 
-1. 一个从头开始重建的 **通用物理引擎**，能够模拟各种材料和物理现象。
-2. 一个 **轻量级**、**超快**、**Python 风格** 和 **用户友好** 的机器人模拟平台。
-3. 一个强大且快速的 **照片级真实感渲染系统**。
-4. 一个 **生成数据引擎**，将用户提示的自然语言描述转换为各种数据模式。
+1. [概述](#概述)
+2. [主要特点](#主要特点)
+3. [快速入门](#快速入门)
+4. [参与贡献](#参与贡献)
+5. [帮助支持](#帮助支持)
+6. [许可证与致谢](#许可证和致谢)
+7. [相关论文](#genesis-背后的论文)
+8. [引用](#引用)
 
-Genesis 由一个重新设计和重建的通用物理引擎驱动，集成了各种物理求解器及其耦合到一个统一的框架中。这个核心物理引擎通过一个生成代理框架在上层进行增强，旨在实现机器人及其他领域的完全 **自动化数据生成**。
+## 概述
 
-目前，我们正在开源 **底层物理引擎和模拟平台**。我们的生成框架是一个模块化系统，包含许多不同的生成模块，每个模块处理一定范围的数据模式，由一个高级代理路由。一些模块集成了现有的论文，一些仍在提交中。我们的生成功能将逐步开放。如果您感兴趣，请参阅下面的 [论文列表](#papers-behind-genesis)。
+Genesis 是专为 *机器人/嵌入式 AI/物理 AI* 应用设计的通用物理平台，集成了以下核心功能：
 
-Genesis 的构建和持续发展基于以下 ***长期使命***：
+- **通用物理引擎**: 从底层重建,支持多种材料和物理现象模拟
+- **机器人模拟平台**: 轻量、高速、Python友好的开发环境
+- **真实感渲染**: 内置光线追踪渲染系统
+- **生成数据引擎**: 自然语言驱动的多模态数据生成
 
-1. **降低使用物理模拟的门槛**，使机器人研究对所有人都可访问。（参见我们的 [承诺](https://genesis-world.readthedocs.io/en/latest/user_guide/overview/mission.html)）
-2. **将各种最先进的物理求解器统一到一个框架中**，使用最先进的模拟技术在虚拟领域中重新创建整个物理世界，达到最高的物理、视觉和感官保真度。
-3. **最小化人类在收集和生成机器人及其他领域数据上的努力**，让数据飞轮自行旋转。
+我们的长期使命:
 
-项目页面：<https://genesis-embodied-ai.github.io/>
+- 降低物理模拟使用门槛
+- 统一各类物理求解器
+- 实现数据生成自动化
+
+项目主页: <https://genesis-embodied-ai.github.io/>
 
 ## 主要特点
 
@@ -43,29 +57,50 @@ Genesis 的构建和持续发展基于以下 ***长期使命***：
 - **基于物理的触觉传感器**：Genesis 包含一个基于物理的可微分 [触觉传感器模拟模块](https://github.com/Genesis-Embodied-AI/DiffTactile)。这将很快集成到公共版本中（预计在 0.3.0 版本中）。
 - **用户友好性**：Genesis 设计为尽可能简化模拟的使用。从安装到 API 设计，如果有任何您觉得不直观或难以使用的地方，请 [告诉我们](https://github.com/Genesis-Embodied-AI/Genesis/issues)。
 
-## 入门
+## 快速入门
 
-### 快速安装
+### 安装
 
 Genesis 可通过 PyPI 获取：
 
 ```bash
-pip install genesis-world  # 需要 Python >=3.9;
+pip install genesis-world  # 需要 Python >=3.9
+```
+
+同时需要按照[官方指南](https://pytorch.org/get-started/locally/)安装 PyTorch。
+
+### Docker 支持
+
+如果您想通过 Docker 使用 Genesis，您可以首先构建 Docker 镜像，命令如下：
+
+```bash
+docker build -t genesis -f docker/Dockerfile docker
 ```
 
-您还需要按照 [官方说明](https://pytorch.org/get-started/locally/) 安装 **PyTorch**。
+然后，您可以在 Docker 镜像内运行示例代码（挂载到 `/workspace/examples`）：
+
+```bash
+xhost +local:root # 允许容器访问显示器
+
+docker run --gpus all --rm -it \
+-e DISPLAY=$DISPLAY \
+-v /tmp/.X11-unix/:/tmp/.X11-unix \
+-v $PWD:/workspace \
+genesis
+```
 
 ### 文档
 
-请参阅我们的 [文档网站（英文）](https://genesis-world.readthedocs.io/en/latest/user_guide/index.html)/[（中文）](https://genesis-world.readthedocs.io/zh-cn/latest/user_guide/index.html)以获取详细的安装步骤、教程和 API 参考。
+- [英文文档](https://genesis-world.readthedocs.io/en/latest/user_guide/index.html)
+- [中文文档](https://genesis-world.readthedocs.io/zh-cn/latest/user_guide/index.html)
 
-## 贡献 Genesis
+## 参与贡献
 
 Genesis 项目的目标是构建一个完全透明、用户友好的生态系统，让来自机器人和计算机图形学的贡献者 **共同创建一个高效、真实（物理和视觉上）的虚拟世界，用于机器人研究及其他领域**。
 
 我们真诚地欢迎来自社区的 *任何形式的贡献*，以使世界对机器人更友好。从 **新功能的拉取请求**、**错误报告**，到甚至是使 Genesis API 更直观的微小 **建议**，我们都全心全意地感谢！
 
-## 支持
+## 帮助支持
 
 - 请使用 Github [Issues](https://github.com/Genesis-Embodied-AI/Genesis/issues) 报告错误和提出功能请求。
 
@@ -116,7 +151,7 @@ Genesis 是一个大规模的努力，将各种现有和正在进行的研究工
 
 如果您在研究中使用了 Genesis，我们将非常感谢您引用它。我们仍在撰写技术报告，在其公开之前，您可以考虑引用：
 
-```
+```bibtex
 @software{Genesis,
   author = {Genesis Authors},
   title = {Genesis: A Universal and Generative Physics Engine for Robotics and Beyond},
diff --git a/README_JA.md b/README_JA.md
new file mode 100644
index 00000000..254ec07c
--- /dev/null
+++ b/README_JA.md
@@ -0,0 +1,171 @@
+![Genesis](imgs/big_text.png)
+
+![Teaser](imgs/teaser.png)
+
+[![PyPI - Version](https://img.shields.io/pypi/v/genesis-world)](https://pypi.org/project/genesis-world/)
+[![PyPI - Downloads](https://img.shields.io/pypi/dm/genesis-world)](https://pypi.org/project/genesis-world/)
+[![GitHub Issues](https://img.shields.io/github/issues/Genesis-Embodied-AI/Genesis)](https://github.com/Genesis-Embodied-AI/Genesis/issues)
+[![GitHub Discussions](https://img.shields.io/github/discussions/Genesis-Embodied-AI/Genesis)](https://github.com/Genesis-Embodied-AI/Genesis/discussions)
+
+[![README in English](https://img.shields.io/badge/English-d9d9d9)](./README.md)
+[![简体中文版自述文件](https://img.shields.io/badge/简体中文-d9d9d9)](./README_CN.md)
+[![日本語版 README](https://img.shields.io/badge/日本語-d9d9d9)](./README_JA.md)
+
+# Genesis
+## 🔥 最新情報
+- [2024-12-25] [レイトレーシングレンダラー](#docker)をサポートするDockerを追加しました。
+- [2024-12-24] [Genesisへの貢献方法](https://github.com/Genesis-Embodied-AI/Genesis/blob/main/CONTRIBUTING.md)に関するガイドラインを追加しました。
+
+## 目次
+
+1. [Genesisとは？](#what-is-genesis)
+2. [主な機能](#key-features)
+3. [簡単インストール](#quick-installation)
+4. [Docker](#docker)
+5. [ドキュメント](#documentation)
+6. [Genesisへの貢献](#contributing-to-genesis)
+7. [サポート](#support)
+8. [ライセンスと謝辞](#license-and-acknowledgments)
+9. [関連論文](#associated-papers)
+10. [引用](#citation)
+
+## Genesisとは？
+
+Genesisは、汎用的な*ロボティクス/身体性を持ったAI*アプリケーション向けに設計された物理シミュレーションプラットフォームです。このプラットフォームは以下のような特徴があります：
+
+1. あらゆる種類の材料や物理現象をシミュレート可能な**汎用物理エンジン**。
+2. **軽量**、**超高速**、**Python的**、そして**ユーザーフレンドリー**なロボティクスシミュレーションプラットフォーム。
+3. 高速で強力な**フォトリアリスティックなレンダリングシステム**。
+4. ユーザーの自然言語による指示をもとに様々なデータモダリティを生成する**生成型データエンジン**。
+
+Genesisの目指すところ：
+
+- **物理シミュレーションのハードルを下げ**、ロボティクス研究を誰でもアクセス可能にすること。詳細は[ミッションステートメント](https://genesis-world.readthedocs.io/en/latest/user_guide/overview/mission.html)をご覧ください。
+- **多様な物理ソルバーを統合**し、最高の忠実度で物理世界を再現すること。
+- **データ生成を自動化**し、人間の労力を削減し、データ生成の効率を最大化すること。
+
+プロジェクトページ: <https://genesis-embodied-ai.github.io/>
+
+## 主な機能
+
+- **速度**: RTX 4090単体でフランカロボットアームを4300万FPS（リアルタイムの43万倍速）でシミュレーション可能。
+- **クロスプラットフォーム**: Linux、macOS、Windowsで動作し、CPU、Nvidia/AMD GPU、Apple Metalをサポート。
+- **多様な物理ソルバーの統合**: 剛体、MPM、SPH、FEM、PBD、安定流体シミュレーション。
+- **幅広い材料モデル**: 剛体、液体、気体、変形体、薄膜オブジェクト、粒状材料などをシミュレーション可能。
+- **様々なロボットへの対応**: ロボットアーム、脚付きロボット、ドローン、*ソフトロボット*など。また、`MJCF (.xml)`、`URDF`、`.obj`、`.glb`、`.ply`、`.stl`などの形式をサポート。
+- **フォトリアルなレンダリング**: レイトレーシングベースのレンダリングをネイティブでサポート。
+- **微分可能性**: 完全な微分可能性を備えた設計。現時点では、MPMソルバーとツールソルバーが対応しており、将来的には他のソルバーも対応予定（まず剛体および連結体ソルバーから開始）。
+- **物理ベースの触覚シミュレーション**: 微分可能な[触覚センサーシミュレーション](https://github.com/Genesis-Embodied-AI/DiffTactile)が近日公開予定（バージョン0.3.0を予定）。
+- **ユーザーフレンドリー**: シンプルで直感的なインストールとAPI設計。
+
+## インストール
+
+GenesisはPyPIで利用可能です：
+
+```bash
+pip install genesis-world  # Python >=3.9 が必要です;
+```
+
+また、**PyTorch**を[公式手順](https://pytorch.org/get-started/locally/)に従ってインストールする必要があります。
+
+最新バージョンを利用するには、リポジトリをクローンしてローカルにインストールしてください：
+
+```bash
+git clone https://github.com/Genesis-Embodied-AI/Genesis.git
+cd Genesis
+pip install -e .
+```
+
+## Docker
+
+DockerからGenesisを利用する場合は、まずDockerイメージをビルドします：
+
+```bash
+docker build -t genesis -f docker/Dockerfile docker
+```
+
+その後、Dockerイメージ内で例を実行できます（`/workspace/examples`にマウント）：
+
+```bash
+xhost +local:root # コンテナがディスプレイにアクセスできるようにする
+
+docker run --gpus all --rm -it \
+-e DISPLAY=$DISPLAY \
+-v /tmp/.X11-unix/:/tmp/.X11-unix \
+-v $PWD:/workspace \
+genesis
+```
+
+## ドキュメント
+
+包括的なドキュメントは現時点では[英語](https://genesis-world.readthedocs.io/en/latest/user_guide/index.html)および[中国語](https://genesis-world.readthedocs.io/zh-cn/latest/user_guide/index.html)で提供されています。詳細なインストール手順、チュートリアル、APIリファレンスが含まれています。
+
+## Genesisへの貢献
+
+Genesisプロジェクトはオープンで協力的な取り組みです。以下を含む、コミュニティからのあらゆる貢献を歓迎します：
+
+- 新機能やバグ修正のための**プルリクエスト**。
+- GitHub Issuesを通じた**バグ報告**。
+- Genesisの使いやすさを向上させるための**提案**。
+
+詳細は[貢献ガイド](https://github.com/Genesis-Embodied-AI/Genesis/blob/main/CONTRIBUTING.md)をご参照ください。
+
+## サポート
+
+- バグ報告や機能リクエストはGitHubの[Issues](https://github.com/Genesis-Embodied-AI/Genesis/issues)をご利用ください。
+- 議論や質問はGitHubの[Discussions](https://github.com/Genesis-Embodied-AI/Genesis/discussions)で行えます。
+
+## ライセンスと謝辞
+
+GenesisのソースコードはApache 2.0ライセンスで提供されています。
+
+Genesisの開発は以下のオープンソースプロジェクトのおかげで可能になりました：
+
+- [Taichi](https://github.com/taichi-dev/taichi): 高性能でクロスプラットフォーム対応の計算バックエンド。Taichiチームの技術サポートに感謝します！
+- [FluidLab](https://github.com/zhouxian/FluidLab): 参照用のMPMソルバー実装。  
+- [SPH_Taichi](https://github.com/erizmr/SPH_Taichi): 参照用のSPHソルバー実装。  
+- [Ten Minute Physics](https://matthias-research.github.io/pages/tenMinutePhysics/index.html) と [PBF3D](https://github.com/WASD4959/PBF3D): 参照用のPBD（粒子ベースの物理）ソルバー実装。  
+- [MuJoCo](https://github.com/google-deepmind/mujoco): 剛体ダイナミクスの参照用実装。  
+- [libccd](https://github.com/danfis/libccd): 衝突検出の参照用実装。  
+- [PyRender](https://github.com/mmatl/pyrender): ラスタライズベースのレンダラー。  
+- [LuisaCompute](https://github.com/LuisaGroup/LuisaCompute) と [LuisaRender](https://github.com/LuisaGroup/LuisaRender): レイトレーシングDSL。
+
+## 関連論文
+
+Genesisプロジェクトに関与した主要な研究論文の一覧：
+
+- Xian, Zhou, et al. "Fluidlab: A differentiable environment for benchmarking complex fluid manipulation." arXiv preprint arXiv:2303.02346 (2023).
+- Xu, Zhenjia, et al. "Roboninja: Learning an adaptive cutting policy for multi-material objects." arXiv preprint arXiv:2302.11553 (2023).
+- Wang, Yufei, et al. "Robogen: Towards unleashing infinite data for automated robot learning via generative simulation." arXiv preprint arXiv:2311.01455 (2023).
+- Wang, Tsun-Hsuan, et al. "Softzoo: A soft robot co-design benchmark for locomotion in diverse environments." arXiv preprint arXiv:2303.09555 (2023).
+- Wang, Tsun-Hsuan Johnson, et al. "Diffusebot: Breeding soft robots with physics-augmented generative diffusion models." Advances in Neural Information Processing Systems 36 (2023): 44398-44423.
+- Katara, Pushkal, Zhou Xian, and Katerina Fragkiadaki. "Gen2sim: Scaling up robot learning in simulation with generative models." 2024 IEEE International Conference on Robotics and Automation (ICRA). IEEE, 2024.
+- Si, Zilin, et al. "DiffTactile: A Physics-based Differentiable Tactile Simulator for Contact-rich Robotic Manipulation." arXiv preprint arXiv:2403.08716 (2024).
+- Wang, Yian, et al. "Thin-Shell Object Manipulations With Differentiable Physics Simulations." arXiv preprint arXiv:2404.00451 (2024).
+- Lin, Chunru, et al. "UBSoft: A Simulation Platform for Robotic Skill Learning in Unbounded Soft Environments." arXiv preprint arXiv:2411.12711 (2024).
+- Zhou, Wenyang, et al. "EMDM: Efficient motion diffusion model for fast and high-quality motion generation." European Conference on Computer Vision. Springer, Cham, 2025.
+- Qiao, Yi-Ling, Junbang Liang, Vladlen Koltun, and Ming C. Lin. "Scalable differentiable physics for learning and control." International Conference on Machine Learning. PMLR, 2020.
+- Qiao, Yi-Ling, Junbang Liang, Vladlen Koltun, and Ming C. Lin. "Efficient differentiable simulation of articulated bodies." In International Conference on Machine Learning, PMLR, 2021.
+- Qiao, Yi-Ling, Junbang Liang, Vladlen Koltun, and Ming Lin. "Differentiable simulation of soft multi-body systems." Advances in Neural Information Processing Systems 34 (2021).
+- Wan, Weilin, et al. "Tlcontrol: Trajectory and language control for human motion synthesis." arXiv preprint arXiv:2311.17135 (2023).
+- Wang, Yian, et al. "Architect: Generating Vivid and Interactive 3D Scenes with Hierarchical 2D Inpainting." arXiv preprint arXiv:2411.09823 (2024).
+- Zheng, Shaokun, et al. "LuisaRender: A high-performance rendering framework with layered and unified interfaces on stream architectures." ACM Transactions on Graphics (TOG) 41.6 (2022): 1-19.
+- Fan, Yingruo, et al. "Faceformer: Speech-driven 3d facial animation with transformers." Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. 2022.
+- Wu, Sichun, Kazi Injamamul Haque, and Zerrin Yumak. "ProbTalk3D: Non-Deterministic Emotion Controllable Speech-Driven 3D Facial Animation Synthesis Using VQ-VAE." Proceedings of the 17th ACM SIGGRAPH Conference on Motion, Interaction, and Games. 2024.
+- Dou, Zhiyang, et al. "C· ase: Learning conditional adversarial skill embeddings for physics-based characters." SIGGRAPH Asia 2023 Conference Papers. 2023.
+
+さらに多数の現在進行形のプロジェクトがあります。
+
+## 引用
+
+研究でGenesisを使用する場合、以下を引用してください：
+
+```bibtex
+@software{Genesis,
+  author = {Genesis Authors},
+  title = {Genesis: A Universal and Generative Physics Engine for Robotics and Beyond},
+  month = {December},
+  year = {2024},
+  url = {https://github.com/Genesis-Embodied-AI/Genesis}
+}
+```
\ No newline at end of file
diff --git a/docker/10_nvidia.json b/docker/10_nvidia.json
new file mode 100644
index 00000000..2bfcca05
--- /dev/null
+++ b/docker/10_nvidia.json
@@ -0,0 +1,6 @@
+{
+    "file_format_version" : "1.0.0",
+    "ICD" : {
+        "library_path" : "libEGL_nvidia.so.0"
+    }
+}
diff --git a/docker/Dockerfile b/docker/Dockerfile
new file mode 100644
index 00000000..e2ef6707
--- /dev/null
+++ b/docker/Dockerfile
@@ -0,0 +1,118 @@
+ARG CUDA_VERSION=12.1
+
+# ===============================================================
+# Stage 1: Build LuisaRender
+# ===============================================================
+FROM pytorch/pytorch:2.5.1-cuda${CUDA_VERSION}-cudnn9-devel AS builder
+
+ENV DEBIAN_FRONTEND=noninteractive
+ARG PYTHON_VERSION=3.11
+
+# Install necessary packages
+RUN apt-get update && apt-get install -y --no-install-recommends \
+    build-essential \
+    manpages-dev \
+    libvulkan-dev \
+    zlib1g-dev \
+    xorg-dev libglu1-mesa-dev \
+    libsnappy-dev \
+    software-properties-common \
+    git \
+    curl \
+    wget
+RUN add-apt-repository ppa:ubuntu-toolchain-r/test && \
+    apt update && \
+    apt install -y --no-install-recommends \
+    gcc-11 \
+    g++-11 \
+    gcc-11 g++-11 patchelf && \
+    rm -rf /var/lib/apt/lists/*
+
+# Set GCC-11 and G++-11 as the default
+RUN update-alternatives --install /usr/bin/gcc gcc /usr/bin/gcc-11 110 && \
+    update-alternatives --install /usr/bin/g++ g++ /usr/bin/g++-11 110
+
+# Install Rust for build requirements
+RUN curl --proto '=https' --tlsv1.2 -sSf https://sh.rustup.rs | sh -s -- -y
+
+RUN pip install "pybind11[global]"
+
+# Install CMake
+RUN wget https://github.com/Kitware/CMake/releases/download/v3.31.0-rc2/cmake-3.31.0-rc2-linux-x86_64.sh && \
+    chmod +x cmake-3.31.0-rc2-linux-x86_64.sh && \
+    ./cmake-3.31.0-rc2-linux-x86_64.sh --skip-license --prefix=/usr/local && \
+    rm cmake-3.31.0-rc2-linux-x86_64.sh
+
+# Build LuisaRender
+WORKDIR /workspace
+RUN git clone https://github.com/Genesis-Embodied-AI/Genesis.git && \
+    cd Genesis && \
+    git submodule update --init --recursive
+COPY build_luisa.sh /workspace/build_luisa.sh
+RUN chmod +x ./build_luisa.sh && ./build_luisa.sh ${PYTHON_VERSION}
+
+# ===============================================================
+# Stage 2: Runtime Environment
+# ===============================================================
+FROM pytorch/pytorch:2.5.1-cuda${CUDA_VERSION}-cudnn9-devel
+
+ARG PYTHON_VERSION=3.11
+ENV DEBIAN_FRONTEND=noninteractive
+ENV NVIDIA_DRIVER_CAPABILITIES=all
+
+# Install runtime dependencies
+RUN apt-get update && apt-get install -y --no-install-recommends \
+    tmux \
+    git \
+    curl \
+    wget \
+    bash-completion \
+    libgl1 \
+    libgl1-mesa-glx \
+    libegl-dev \
+    libegl1 \
+    libxrender1 \
+    libglib2.0-0 \
+    ffmpeg \
+    libgtk2.0-dev \
+    pkg-config \
+    libvulkan-dev \
+    libgles2 \
+    libglvnd0 \
+    libglx0 \
+    && apt clean \
+    && rm -rf /var/lib/apt/lists/*
+
+WORKDIR /workspace
+
+# --------------------------- Genesis ----------------------------
+RUN pip install --no-cache-dir open3d
+RUN git clone https://github.com/Genesis-Embodied-AI/Genesis.git && \
+    cd Genesis && \
+    pip install . && \
+    pip install --no-cache-dir PyOpenGL==3.1.5
+
+# ------------------------ Motion planning -----------------------
+RUN PYTHON_MAJOR_MINOR=$(echo ${PYTHON_VERSION} | tr -d '.') && \
+    wget https://github.com/ompl/ompl/releases/download/prerelease/ompl-1.6.0-cp${PYTHON_MAJOR_MINOR}-cp${PYTHON_MAJOR_MINOR}-manylinux_2_28_x86_64.whl && \
+    pip install ompl-1.6.0-cp${PYTHON_MAJOR_MINOR}-cp${PYTHON_MAJOR_MINOR}-manylinux_2_28_x86_64.whl && \
+    rm ompl-1.6.0-cp${PYTHON_MAJOR_MINOR}-cp${PYTHON_MAJOR_MINOR}-manylinux_2_28_x86_64.whl
+
+# -------------------- Surface Reconstruction --------------------
+# Set the LD_LIBRARY_PATH directly in the environment
+COPY --from=builder /workspace/Genesis/genesis/ext/ParticleMesher/ParticleMesherPy /opt/conda/lib/python3.1/site-packages/genesis/ext/ParticleMesher/ParticleMesherPy
+ENV LD_LIBRARY_PATH=/opt/conda/lib/python3.1/site-packages/genesis/ext/ParticleMesher/ParticleMesherPy:$LD_LIBRARY_PATH
+
+# --------------------- Ray Tracing Renderer ---------------------
+# Copy LuisaRender build artifacts from the builder stage
+COPY --from=builder /workspace/Genesis/genesis/ext/LuisaRender/build/bin /opt/conda/lib/python3.1/site-packages/genesis/ext/LuisaRender/build/bin
+# fix GLIBCXX_3.4.30 not found
+RUN cd /opt/conda/lib && \
+    mv libstdc++.so.6 libstdc++.so.6.old && \
+    ln -s /usr/lib/x86_64-linux-gnu/libstdc++.so.6 libstdc++.so.6
+
+COPY 10_nvidia.json /usr/share/glvnd/egl_vendor.d/10_nvidia.json
+COPY nvidia_icd.json /usr/share/vulkan/icd.d/nvidia_icd.json
+COPY nvidia_layers.json /etc/vulkan/implicit_layer.d/nvidia_layers.json
+
+ENTRYPOINT ["/bin/bash"]
diff --git a/docker/build_luisa.sh b/docker/build_luisa.sh
new file mode 100644
index 00000000..95d861c5
--- /dev/null
+++ b/docker/build_luisa.sh
@@ -0,0 +1,22 @@
+#!/bin/bash
+
+# Check if Python version is provided
+if [ -z "$1" ]; then
+    echo "Usage: $0 <PYTHON_VERSION>"
+    exit 1
+fi
+
+PYTHON_VERSION=$1
+
+cd Genesis/genesis/ext/LuisaRender && \
+git submodule update --init --recursive && \
+mkdir -p build && \
+cmake -S . -B build \
+    -D CMAKE_BUILD_TYPE=Release \
+    -D PYTHON_VERSIONS=$PYTHON_VERSION \
+    -D LUISA_COMPUTE_DOWNLOAD_NVCOMP=ON \
+    -D LUISA_COMPUTE_DOWNLOAD_OIDN=ON \
+    -D LUISA_COMPUTE_ENABLE_GUI=OFF \
+    -D LUISA_COMPUTE_ENABLE_CUDA=ON \
+    -Dpybind11_DIR=$(python3 -c "import pybind11; print(pybind11.get_cmake_dir())") && \
+cmake --build build -j $(nproc)
\ No newline at end of file
diff --git a/docker/nvidia_icd.json b/docker/nvidia_icd.json
new file mode 100644
index 00000000..69600b17
--- /dev/null
+++ b/docker/nvidia_icd.json
@@ -0,0 +1,7 @@
+{
+    "file_format_version" : "1.0.0",
+    "ICD": {
+        "library_path": "libGLX_nvidia.so.0",
+        "api_version" : "1.2.155"
+    }
+}
diff --git a/docker/nvidia_layers.json b/docker/nvidia_layers.json
new file mode 100644
index 00000000..a8e098eb
--- /dev/null
+++ b/docker/nvidia_layers.json
@@ -0,0 +1,22 @@
+
+{
+    "file_format_version" : "1.0.0",
+    "layer": {
+        "name": "VK_LAYER_NV_optimus",
+        "type": "INSTANCE",
+        "library_path": "libGLX_nvidia.so.0",
+        "api_version" : "1.2.155",
+        "implementation_version" : "1",
+        "description" : "NVIDIA Optimus layer",
+        "functions": {
+            "vkGetInstanceProcAddr": "vk_optimusGetInstanceProcAddr",
+            "vkGetDeviceProcAddr": "vk_optimusGetDeviceProcAddr"
+        },
+        "enable_environment": {
+            "__NV_PRIME_RENDER_OFFLOAD": "1"
+        },
+        "disable_environment": {
+            "DISABLE_LAYER_NV_OPTIMUS_1": ""
+        }
+    }
+}
diff --git a/examples/drone/README.md b/examples/drone/README.md
index 349a8d18..39eb37c1 100644
--- a/examples/drone/README.md
+++ b/examples/drone/README.md
@@ -27,6 +27,50 @@ Run with:
 python fly.py -v -m
 ```
 
+### 3. Hover Environment (`hover_env.py`, `hover_train.py`, `hover_eval.py`)
+
+The hover environment (`hover_env.py`) is designed to train a drone to maintain a stable hover position by reaching randomly generated target points. The environment includes:
+
+ - Initialization of the scene and entities (plane, drone and target).
+ - Reward functions to provide feedback to the agent based on its performance in reaching the target points.
+ - **Command resampling to generate new random target points** and environment reset functionalities to ensure continuous training.
+
+**Acknowledgement**: The reward design is inspired by [Champion-level drone racing using deep
+reinforcement learning (Nature 2023)](https://www.nature.com/articles/s41586-023-06419-4.pdf)
+
+#### 3.0 Installation
+
+At this stage, we have defined the environments. Now, we use the PPO implementation from `rsl-rl` to train the policy. Follow these installation steps:
+
+```bash
+# Install rsl_rl.
+git clone https://github.com/leggedrobotics/rsl_rl
+cd rsl_rl && git checkout v1.0.2 && pip install -e .
+
+# Install tensorboard.
+pip install tensorboard
+```
+
+#### 3.1 Training
+
+Train the drone hovering policy using the `HoverEnv` environment.
+
+Run with:
+
+```bash
+python hover_train.py -e drone-hovering -B 8192 --max_iterations 500
+```
+
+#### 3.2 Evaluation
+
+Evaluate the trained drone hovering policy.
+
+Run with:
+
+```bash
+python hover_eval.py -e drone-hovering --ckpt 500 --record
+```
+
 ## Technical Details
 
 - The drone model used is the Crazyflie 2.X (`urdf/drones/cf2x.urdf`)
diff --git a/examples/drone/hover_env.py b/examples/drone/hover_env.py
new file mode 100644
index 00000000..3d48d6b3
--- /dev/null
+++ b/examples/drone/hover_env.py
@@ -0,0 +1,260 @@
+import torch
+import math
+import genesis as gs
+from genesis.utils.geom import quat_to_xyz, transform_by_quat, inv_quat, transform_quat_by_quat
+
+def gs_rand_float(lower, upper, shape, device):
+    return (upper - lower) * torch.rand(size=shape, device=device) + lower
+
+class HoverEnv:
+    def __init__(self, num_envs, env_cfg, obs_cfg, reward_cfg, command_cfg, show_viewer=False, device="cuda"):
+        self.device = torch.device(device)
+
+        self.num_envs = num_envs
+        self.num_obs = obs_cfg["num_obs"]
+        self.num_privileged_obs = None
+        self.num_actions = env_cfg["num_actions"]
+        self.num_commands = command_cfg["num_commands"]
+
+        self.simulate_action_latency = env_cfg["simulate_action_latency"]
+        self.dt = 0.01  # run in 100hz
+        self.max_episode_length = math.ceil(env_cfg["episode_length_s"] / self.dt)
+
+        self.env_cfg = env_cfg
+        self.obs_cfg = obs_cfg
+        self.reward_cfg = reward_cfg
+        self.command_cfg = command_cfg
+
+        self.obs_scales = obs_cfg["obs_scales"]
+        self.reward_scales = reward_cfg["reward_scales"]
+
+        # create scene
+        self.scene = gs.Scene(
+            sim_options=gs.options.SimOptions(dt=self.dt, substeps=2),
+            viewer_options=gs.options.ViewerOptions(
+                max_FPS=env_cfg["max_visualize_FPS"],
+                camera_pos=(3.0, 0.0, 3.0),
+                camera_lookat=(0.0, 0.0, 1.0),
+                camera_fov=40,
+            ),
+            vis_options=gs.options.VisOptions(n_rendered_envs=1),
+            rigid_options=gs.options.RigidOptions(
+                dt=self.dt,
+                constraint_solver=gs.constraint_solver.Newton,
+                enable_collision=True,
+                enable_joint_limit=True,
+            ),
+            show_viewer=show_viewer,
+        )
+
+        # add plane
+        self.scene.add_entity(gs.morphs.Plane())
+
+        # add target
+        if self.env_cfg["visualize_target"]:
+            self.target = self.scene.add_entity(morph=gs.morphs.Mesh(
+                                    file="meshes/sphere.obj",
+                                    scale=0.05, 
+                                    fixed=True, 
+                                    collision=False,
+                                ),
+                                surface=gs.surfaces.Rough(
+                                    diffuse_texture=gs.textures.ColorTexture(
+                                        color=(1.0, 0.5, 0.5),
+                                    ),
+                                ),
+                            )
+
+        # add camera
+        if self.env_cfg["visualize_camera"]:
+            self.cam = self.scene.add_camera(
+                res=(640, 480),
+                pos=(3.5, 0.0, 2.5),
+                lookat=(0, 0, 0.5),
+                fov=30,
+                GUI=True,
+            )
+
+        # add drone
+        self.base_init_pos = torch.tensor(self.env_cfg["base_init_pos"], device=self.device)
+        self.base_init_quat = torch.tensor(self.env_cfg["base_init_quat"], device=self.device)
+        self.inv_base_init_quat = inv_quat(self.base_init_quat)
+        self.drone = self.scene.add_entity(gs.morphs.Drone(file="urdf/drones/cf2x.urdf"))
+
+        # build scene
+        self.scene.build(n_envs=num_envs)
+
+        # prepare reward functions and multiply reward scales by dt
+        self.reward_functions, self.episode_sums = dict(), dict()
+        for name in self.reward_scales.keys():
+            self.reward_scales[name] *= self.dt
+            self.reward_functions[name] = getattr(self, "_reward_" + name)
+            self.episode_sums[name] = torch.zeros((self.num_envs,), device=self.device, dtype=gs.tc_float)
+
+        # initialize buffers
+        self.obs_buf = torch.zeros((self.num_envs, self.num_obs), device=self.device, dtype=gs.tc_float)
+        self.rew_buf = torch.zeros((self.num_envs,), device=self.device, dtype=gs.tc_float)
+        self.reset_buf = torch.ones((self.num_envs,), device=self.device, dtype=gs.tc_int)
+        self.episode_length_buf = torch.zeros((self.num_envs,), device=self.device, dtype=gs.tc_int)
+        self.commands = torch.zeros((self.num_envs, self.num_commands), device=self.device, dtype=gs.tc_float)
+
+        self.actions = torch.zeros((self.num_envs, self.num_actions), device=self.device, dtype=gs.tc_float)
+        self.last_actions = torch.zeros_like(self.actions)
+
+        self.base_pos = torch.zeros((self.num_envs, 3), device=self.device, dtype=gs.tc_float)
+        self.base_quat = torch.zeros((self.num_envs, 4), device=self.device, dtype=gs.tc_float)
+        self.base_lin_vel = torch.zeros((self.num_envs, 3), device=self.device, dtype=gs.tc_float)
+        self.base_ang_vel = torch.zeros((self.num_envs, 3), device=self.device, dtype=gs.tc_float)
+        self.last_base_pos = torch.zeros_like(self.base_pos)
+
+        self.extras = dict()  # extra information for logging
+
+    def _resample_commands(self, envs_idx):
+        self.commands[envs_idx, 0] = gs_rand_float(*self.command_cfg["pos_x_range"], (len(envs_idx),), self.device)
+        self.commands[envs_idx, 1] = gs_rand_float(*self.command_cfg["pos_y_range"], (len(envs_idx),), self.device)
+        self.commands[envs_idx, 2] = gs_rand_float(*self.command_cfg["pos_z_range"], (len(envs_idx),), self.device)
+        if self.target is not None:
+            self.target.set_pos(self.commands[envs_idx], zero_velocity=True, envs_idx=envs_idx)
+
+    def _at_target(self):
+        at_target = (
+            (torch.norm(self.rel_pos, dim=1) < self.env_cfg["at_target_threshold"])
+            .nonzero(as_tuple=False)
+            .flatten()
+        )
+        return at_target
+
+    def step(self, actions):
+        self.actions = torch.clip(actions, -self.env_cfg["clip_actions"], self.env_cfg["clip_actions"])
+        exec_actions = self.actions.cpu()
+        # exec_actions = self.last_actions.cpu() if self.simulate_action_latency else self.actions.cpu()
+        # target_dof_pos = exec_actions * self.env_cfg["action_scale"] + self.default_dof_pos
+        # self.drone.control_dofs_position(target_dof_pos)
+
+        # 14468 is hover rpm
+        self.drone.set_propellels_rpm((1 + exec_actions*0.8) * 14468.429183500699)
+        self.scene.step()
+
+        # update buffers
+        self.episode_length_buf += 1
+        self.last_base_pos[:] = self.base_pos[:]
+        self.base_pos[:] = self.drone.get_pos()
+        self.rel_pos = self.commands - self.base_pos
+        self.last_rel_pos = self.commands - self.last_base_pos
+        self.base_quat[:] = self.drone.get_quat()
+        self.base_euler = quat_to_xyz(
+            transform_quat_by_quat(torch.ones_like(self.base_quat) * self.inv_base_init_quat, self.base_quat)
+        )
+        inv_base_quat = inv_quat(self.base_quat)
+        self.base_lin_vel[:] = transform_by_quat(self.drone.get_vel(), inv_base_quat)
+        self.base_ang_vel[:] = transform_by_quat(self.drone.get_ang(), inv_base_quat)
+
+        # resample commands
+        envs_idx = self._at_target()
+        self._resample_commands(envs_idx)
+
+        # check termination and reset
+        self.crash_condition = (
+            (torch.abs(self.base_euler[:, 1]) > self.env_cfg["termination_if_pitch_greater_than"]) |
+            (torch.abs(self.base_euler[:, 0]) > self.env_cfg["termination_if_roll_greater_than"]) |
+            (torch.abs(self.rel_pos[:, 0]) > self.env_cfg["termination_if_x_greater_than"]) |
+            (torch.abs(self.rel_pos[:, 1]) > self.env_cfg["termination_if_y_greater_than"]) |
+            (torch.abs(self.rel_pos[:, 2]) > self.env_cfg["termination_if_z_greater_than"]) |
+            (self.base_pos[:, 2] < self.env_cfg["termination_if_close_to_ground"])
+        )
+        self.reset_buf = (self.episode_length_buf > self.max_episode_length) | self.crash_condition
+
+        time_out_idx = (self.episode_length_buf > self.max_episode_length).nonzero(as_tuple=False).flatten()
+        self.extras["time_outs"] = torch.zeros_like(self.reset_buf, device=self.device, dtype=gs.tc_float)
+        self.extras["time_outs"][time_out_idx] = 1.0
+
+        self.reset_idx(self.reset_buf.nonzero(as_tuple=False).flatten())
+
+        # compute reward
+        self.rew_buf[:] = 0.0
+        for name, reward_func in self.reward_functions.items():
+            rew = reward_func() * self.reward_scales[name]
+            self.rew_buf += rew
+            self.episode_sums[name] += rew
+
+        # compute observations
+        self.obs_buf = torch.cat(
+            [
+                torch.clip(self.rel_pos * self.obs_scales["rel_pos"], -1, 1),
+                self.base_quat,
+                torch.clip(self.base_lin_vel * self.obs_scales["lin_vel"], -1, 1),
+                torch.clip(self.base_ang_vel * self.obs_scales["ang_vel"], -1, 1),
+                self.last_actions,
+            ],
+            axis=-1,
+        )
+
+        self.last_actions[:] = self.actions[:]
+
+        return self.obs_buf, None, self.rew_buf, self.reset_buf, self.extras
+
+    def get_observations(self):
+        return self.obs_buf
+
+    def get_privileged_observations(self):
+        return None
+
+    def reset_idx(self, envs_idx):
+        if len(envs_idx) == 0:
+            return
+
+        # reset base
+        self.base_pos[envs_idx] = self.base_init_pos
+        self.last_base_pos[envs_idx] = self.base_init_pos
+        self.rel_pos = self.commands - self.base_pos
+        self.last_rel_pos = self.commands - self.last_base_pos
+        self.base_quat[envs_idx] = self.base_init_quat.reshape(1, -1)
+        self.drone.set_pos(self.base_pos[envs_idx], zero_velocity=True, envs_idx=envs_idx)
+        self.drone.set_quat(self.base_quat[envs_idx], zero_velocity=True, envs_idx=envs_idx)
+        self.base_lin_vel[envs_idx] = 0
+        self.base_ang_vel[envs_idx] = 0
+        self.drone.zero_all_dofs_velocity(envs_idx)
+
+        # reset buffers
+        self.last_actions[envs_idx] = 0.0
+        self.episode_length_buf[envs_idx] = 0
+        self.reset_buf[envs_idx] = True
+
+        # fill extras
+        self.extras["episode"] = {}
+        for key in self.episode_sums.keys():
+            self.extras["episode"]["rew_" + key] = (
+                torch.mean(self.episode_sums[key][envs_idx]).item() / self.env_cfg["episode_length_s"]
+            )
+            self.episode_sums[key][envs_idx] = 0.0
+
+        self._resample_commands(envs_idx)
+
+    def reset(self):
+        self.reset_buf[:] = True
+        self.reset_idx(torch.arange(self.num_envs, device=self.device))
+        return self.obs_buf, None
+
+    # ------------ reward functions----------------
+    def _reward_target(self):
+        target_rew = torch.sum(torch.square(self.last_rel_pos), dim=1) - torch.sum(torch.square(self.rel_pos), dim=1)
+        return target_rew
+
+    def _reward_smooth(self):
+        smooth_rew = torch.sum(torch.square(self.actions - self.last_actions), dim=1)
+        return smooth_rew
+
+    def _reward_yaw(self):
+        yaw = self.base_euler[:, 2]
+        yaw = torch.where(yaw > 180, yaw - 360, yaw)/180*3.14159    # use rad for yaw_reward
+        yaw_rew = torch.exp(self.reward_cfg["yaw_lambda"] * torch.abs(yaw))
+        return yaw_rew
+    
+    def _reward_angular(self):
+        angular_rew = torch.norm(self.base_ang_vel/3.14159, dim=1)
+        return angular_rew
+
+    def _reward_crash(self):
+        crash_rew = torch.zeros((self.num_envs,), device=self.device, dtype=gs.tc_float)
+        crash_rew[self.crash_condition] = 1
+        return crash_rew
\ No newline at end of file
diff --git a/examples/drone/hover_eval.py b/examples/drone/hover_eval.py
new file mode 100644
index 00000000..9ef11e98
--- /dev/null
+++ b/examples/drone/hover_eval.py
@@ -0,0 +1,68 @@
+import argparse
+import os
+import pickle
+
+import torch
+from hover_env import HoverEnv
+from rsl_rl.runners import OnPolicyRunner
+
+import genesis as gs
+
+
+def main():
+    parser = argparse.ArgumentParser()
+    parser.add_argument("-e", "--exp_name", type=str, default="drone-hovering")
+    parser.add_argument("--ckpt", type=int, default=500)
+    parser.add_argument("--record", action="store_true", default=False)
+    args = parser.parse_args()
+
+    gs.init()
+
+    log_dir = f"logs/{args.exp_name}"
+    env_cfg, obs_cfg, reward_cfg, command_cfg, train_cfg = pickle.load(open(f"logs/{args.exp_name}/cfgs.pkl", "rb"))
+    reward_cfg["reward_scales"] = {}
+
+    # visualize the target
+    env_cfg["visualize_target"] = True
+    # for video recording
+    env_cfg["visualize_camera"] = args.record
+    # set the max FPS for visualization
+    env_cfg["max_visualize_FPS"] = 60
+
+    env = HoverEnv(
+        num_envs=1,
+        env_cfg=env_cfg,
+        obs_cfg=obs_cfg,
+        reward_cfg=reward_cfg,
+        command_cfg=command_cfg,
+        show_viewer=True,
+    )
+
+    runner = OnPolicyRunner(env, train_cfg, log_dir, device="cuda:0")
+    resume_path = os.path.join(log_dir, f"model_{args.ckpt}.pt")
+    runner.load(resume_path)
+    policy = runner.get_inference_policy(device="cuda:0")
+
+    obs, _ = env.reset()
+
+    max_sim_step = int(env_cfg["episode_length_s"]*env_cfg["max_visualize_FPS"])
+    with torch.no_grad():
+        if args.record:
+            env.cam.start_recording()
+            for _ in range(max_sim_step):
+                actions = policy(obs)
+                obs, _, rews, dones, infos = env.step(actions)
+                env.cam.render()
+            env.cam.stop_recording(save_to_filename="video.mp4", fps=env_cfg["max_visualize_FPS"])
+        else:
+            for _ in range(max_sim_step):
+                actions = policy(obs)
+                obs, _, rews, dones, infos = env.step(actions)
+
+if __name__ == "__main__":
+    main()
+
+"""
+# evaluation
+python examples/drone/hover_eval.py
+"""
diff --git a/examples/drone/hover_train.py b/examples/drone/hover_train.py
new file mode 100644
index 00000000..200d42d0
--- /dev/null
+++ b/examples/drone/hover_train.py
@@ -0,0 +1,147 @@
+import argparse
+import os
+import pickle
+import shutil
+
+from hover_env import HoverEnv
+from rsl_rl.runners import OnPolicyRunner
+
+import genesis as gs
+
+
+def get_train_cfg(exp_name, max_iterations):
+
+    train_cfg_dict = {
+        "algorithm": {
+            "clip_param": 0.2,
+            "desired_kl": 0.01,
+            "entropy_coef": 0.002,
+            "gamma": 0.99,
+            "lam": 0.95,
+            "learning_rate": 0.0003,
+            "max_grad_norm": 1.0,
+            "num_learning_epochs": 5,
+            "num_mini_batches": 4,
+            "schedule": "adaptive",
+            "use_clipped_value_loss": True,
+            "value_loss_coef": 1.0,
+        },
+        "init_member_classes": {},
+        "policy": {
+            "activation": "tanh",
+            "actor_hidden_dims": [128, 128],
+            "critic_hidden_dims": [128, 128],
+            "init_noise_std": 1.0,
+        },
+        "runner": {
+            "algorithm_class_name": "PPO",
+            "checkpoint": -1,
+            "experiment_name": exp_name,
+            "load_run": -1,
+            "log_interval": 1,
+            "max_iterations": max_iterations,
+            "num_steps_per_env": 24,
+            "policy_class_name": "ActorCritic",
+            "record_interval": -1,
+            "resume": False,
+            "resume_path": None,
+            "run_name": "",
+            "runner_class_name": "runner_class_name",
+            "save_interval": 100,
+        },
+        "runner_class_name": "OnPolicyRunner",
+        "seed": 1,
+    }
+
+    return train_cfg_dict
+
+
+def get_cfgs():
+    env_cfg = {
+        "num_actions": 4,
+        # termination
+        "termination_if_roll_greater_than": 180,  # degree
+        "termination_if_pitch_greater_than": 180,
+        "termination_if_close_to_ground": 0.1,
+        "termination_if_x_greater_than": 3.0,
+        "termination_if_y_greater_than": 3.0,
+        "termination_if_z_greater_than": 2.0,
+        # base pose
+        "base_init_pos": [0.0, 0.0, 1.0],
+        "base_init_quat": [1.0, 0.0, 0.0, 0.0],
+        "episode_length_s": 15.0,
+        "at_target_threshold": 0.1,
+        "resampling_time_s": 3.0,
+        "simulate_action_latency": True,
+        "clip_actions": 1.0,
+        # visualization
+        "visualize_target": False,
+        "visualize_camera": False,
+        "max_visualize_FPS": 60,
+    }
+    obs_cfg = {
+        "num_obs": 17,
+        "obs_scales": {
+            "rel_pos": 1 / 3.0,
+            "lin_vel": 1 / 3.0,
+            "ang_vel": 1 / 3.14159,
+        },
+    }
+    reward_cfg = {
+        "yaw_lambda": -10.0,
+        "reward_scales": {
+            "target": 10.0,
+            "smooth": -1e-4,
+            "yaw": 0.01,
+            "angular": -2e-4,
+            "crash": -10.0,
+        }
+    }
+    command_cfg = {
+        "num_commands": 3,
+        "pos_x_range": [-1.0, 1.0],
+        "pos_y_range": [-1.0, 1.0],
+        "pos_z_range": [1.0, 1.0],
+    }
+
+    return env_cfg, obs_cfg, reward_cfg, command_cfg
+
+
+def main():
+    parser = argparse.ArgumentParser()
+    parser.add_argument("-e", "--exp_name", type=str, default="drone-hovering")
+    parser.add_argument("-B", "--num_envs", type=int, default=8192)
+    parser.add_argument("--max_iterations", type=int, default=500)
+    args = parser.parse_args()
+
+    gs.init(logging_level="warning")
+
+    log_dir = f"logs/{args.exp_name}"
+    env_cfg, obs_cfg, reward_cfg, command_cfg = get_cfgs()
+    train_cfg = get_train_cfg(args.exp_name, args.max_iterations)
+
+    if os.path.exists(log_dir):
+        shutil.rmtree(log_dir)
+    os.makedirs(log_dir, exist_ok=True)
+
+    env = HoverEnv(
+        num_envs=args.num_envs, env_cfg=env_cfg, obs_cfg=obs_cfg, reward_cfg=reward_cfg, command_cfg=command_cfg
+    )
+
+    runner = OnPolicyRunner(env, train_cfg, log_dir, device="cuda:0")
+
+    pickle.dump(
+        [env_cfg, obs_cfg, reward_cfg, command_cfg, train_cfg],
+        open(f"{log_dir}/cfgs.pkl", "wb"),
+    )
+
+    runner.learn(num_learning_iterations=args.max_iterations, init_at_random_ep_len=True)
+
+
+if __name__ == "__main__":
+    main()
+
+"""
+# training
+python examples/drone/hover_train.py
+"""
diff --git a/examples/rigid/set_phys_attr.py b/examples/rigid/set_phys_attr.py
new file mode 100644
index 00000000..8dc0bf56
--- /dev/null
+++ b/examples/rigid/set_phys_attr.py
@@ -0,0 +1,229 @@
+import argparse
+
+import numpy as np
+
+import genesis as gs
+
+
+def main():
+    parser = argparse.ArgumentParser()
+    parser.add_argument("-v", "--vis", action="store_true", default=False)
+    args = parser.parse_args()
+
+    ########################## init ##########################
+    gs.init(backend=gs.gpu)
+
+    ########################## create a scene ##########################
+    viewer_options = gs.options.ViewerOptions(
+        camera_pos=(0, -3.5, 2.5),
+        camera_lookat=(0.0, 0.0, 0.5),
+        camera_fov=40,
+        max_FPS=60,
+    )
+
+    scene = gs.Scene(
+        viewer_options=viewer_options,
+        sim_options=gs.options.SimOptions(
+            dt=0.01,
+        ),
+        show_viewer=args.vis,
+        rigid_options=gs.options.RigidOptions(
+            # NOTE: Batching dofs/links info to set different physical attributes across environments (in parallel)
+            #       By default, both are False as it's faster and thus only turn this on if necessary
+            batch_dofs_info=True,
+            batch_links_info=True,
+        ),
+    )
+
+    ########################## entities ##########################
+    plane = scene.add_entity(
+        gs.morphs.Plane(),
+    )
+    franka = scene.add_entity(
+        gs.morphs.MJCF(file="xml/franka_emika_panda/panda.xml"),
+    )
+    ########################## build ##########################
+    scene.build(n_envs=2)  # test with 2 different environments
+
+    jnt_names = [
+        "joint1",
+        "joint2",
+        "joint3",
+        "joint4",
+        "joint5",
+        "joint6",
+        "joint7",
+        "finger_joint1",
+        "finger_joint2",
+    ]
+    dofs_idx = [franka.get_joint(name).dof_idx_local for name in jnt_names]
+
+    lnk_names = [
+        "link0",
+        "link1",
+        "link2",
+        "link3",
+        "link4",
+        "link5",
+        "link6",
+        "link7",
+        "hand",
+        "left_finger",
+        "right_finger",
+    ]
+    links_idx = [franka.get_link(name).idx_local for name in lnk_names]
+
+    # Optional: set control gains
+    franka.set_dofs_kp(
+        np.array(
+            [
+                [4500, 4500, 3500, 3500, 2000, 2000, 2000, 100, 100],
+                [100, 100, 2000, 2000, 2000, 3500, 3500, 4500, 4500],
+            ]
+        ),
+        dofs_idx,
+    )
+    print("=== kp ===\n", franka.get_dofs_kp())
+    franka.set_dofs_kv(
+        np.array(
+            [
+                [450, 450, 350, 350, 200, 200, 200, 10, 10],
+                [10, 10, 200, 200, 200, 350, 350, 450, 450],
+            ]
+        ),
+        dofs_idx,
+    )
+    print("=== kv ===\n", franka.get_dofs_kv())
+    franka.set_dofs_force_range(
+        np.array(
+            [
+                [-87, -87, -87, -87, -12, -12, -12, -100, -100],
+                [-120, -100, -12, -12, -12, -87, -87, -87, -87],
+            ]
+        ),
+        np.array(
+            [
+                [87, 87, 87, 87, 12, 12, 12, 100, 100],
+                [100, 100, 12, 12, 12, 87, 87, 87, 87],
+            ]
+        ),
+        dofs_idx,
+    )
+    print("=== force range ===\n", franka.get_dofs_force_range())
+    franka.set_dofs_armature(
+        np.array(
+            [
+                [0.1] * len(dofs_idx),
+                [0.2] * len(dofs_idx),
+            ]
+        ),
+        dofs_idx,
+    )
+    print("=== armature ===\n", franka.get_dofs_armature())
+    franka.set_dofs_stiffness(
+        np.array(
+            [
+                [0.0] * len(dofs_idx),
+                [0.1] * len(dofs_idx),
+            ]
+        ),
+        dofs_idx,
+    )
+    print("=== stiffness ===\n", franka.get_dofs_stiffness())
+    franka.set_dofs_invweight(
+        np.array(
+            [
+                [5.5882, 0.9693, 6.8053, 3.9007, 7.8085, 6.6139, 9.4213, 8.6984, 8.6984],
+                [8.6984, 8.6984, 9.4213, 6.6139, 7.8085, 3.9007, 6.8053, 0.9693, 5.5882],
+            ]
+        ),
+        dofs_idx,
+    )
+    print("=== invweight ===\n", franka.get_dofs_invweight())
+    franka.set_dofs_damping(
+        np.array(
+            [
+                [1.0] * len(dofs_idx),
+                [2.0] * len(dofs_idx),
+            ]
+        ),
+        dofs_idx,
+    )
+    print("=== damping ===\n", franka.get_dofs_damping())
+    franka.set_links_inertial_mass(
+        np.array(
+            [
+                [0.6298, 4.9707, 0.6469, 3.2286, 3.5879, 1.2259, 1.6666, 0.7355, 0.7300, 0.0150, 0.0150],
+                [0.015, 0.015, 0.73, 0.7355, 1.6666, 1.2259, 3.5879, 3.2286, 0.6469, 4.9707, 0.6298],
+            ]
+        ),
+        links_idx,
+    )
+    print("=== links inertial mass ===\n", franka.get_links_inertial_mass())
+    franka.set_links_invweight(
+        np.array(
+            [
+                [0.0, 3.6037e-05, 0.00030664, 0.025365, 0.036351, 0.072328, 0.089559, 0.11661, 0.11288, 3.0179, 3.0179],
+                [3.0179, 3.0179, 0.11288, 0.11661, 0.089559, 0.072328, 0.036351, 0.025365, 0.00030664, 3.6037e-05, 0.0],
+            ]
+        ),
+        links_idx,
+    )
+    print("=== links invweight ===\n", franka.get_links_invweight())
+
+    # Hard reset
+    for i in range(150):
+        if i < 50:
+            franka.set_dofs_position(
+                np.array([1, 1, 0, 0, 0, 0, 0, 0.04, 0.04])[None, :].repeat(scene.n_envs, 0), dofs_idx
+            )
+        elif i < 100:
+            franka.set_dofs_position(
+                np.array([-1, 0.8, 1, -2, 1, 0.5, -0.5, 0.04, 0.04])[None, :].repeat(scene.n_envs, 0), dofs_idx
+            )
+        else:
+            franka.set_dofs_position(np.array([0, 0, 0, 0, 0, 0, 0, 0, 0])[None, :].repeat(scene.n_envs, 0), dofs_idx)
+
+        scene.step()
+
+    # PD control
+    for i in range(1250):
+        if i == 0:
+            franka.control_dofs_position(
+                np.array([1, 1, 0, 0, 0, 0, 0, 0.04, 0.04])[None, :].repeat(scene.n_envs, 0),
+                dofs_idx,
+            )
+        elif i == 250:
+            franka.control_dofs_position(
+                np.array([-1, 0.8, 1, -2, 1, 0.5, -0.5, 0.04, 0.04])[None, :].repeat(scene.n_envs, 0),
+                dofs_idx,
+            )
+        elif i == 500:
+            franka.control_dofs_position(
+                np.array([0, 0, 0, 0, 0, 0, 0, 0, 0])[None, :].repeat(scene.n_envs, 0),
+                dofs_idx,
+            )
+        elif i == 750:
+            # control first dof with velocity, and the rest with position
+            franka.control_dofs_position(
+                np.array([0, 0, 0, 0, 0, 0, 0, 0, 0])[1:][None, :].repeat(scene.n_envs, 0),
+                dofs_idx[1:],
+            )
+            franka.control_dofs_velocity(
+                np.array([1.0, 0, 0, 0, 0, 0, 0, 0, 0])[:1][None, :].repeat(scene.n_envs, 0),
+                dofs_idx[:1],
+            )
+        elif i == 1000:
+            franka.control_dofs_force(
+                np.array([0, 0, 0, 0, 0, 0, 0, 0, 0])[None, :].repeat(scene.n_envs, 0),
+                dofs_idx,
+            )
+        # This is the internal control force computed based on the given control command
+        # If using force control, it's the same as the given control command
+        print("control force:", franka.get_dofs_control_force(dofs_idx))
+
+        scene.step()
+
+
+if __name__ == "__main__":
+    main()
diff --git a/genesis/datatypes.py b/genesis/datatypes.py
index 83960362..c15b0442 100644
--- a/genesis/datatypes.py
+++ b/genesis/datatypes.py
@@ -72,10 +72,6 @@ def _repr_brief(self):
 
         return repr_str
 
-    def __repr__(self):
-        if not __debug__:
-            self.__colorized__repr__()
-
     def __colorized__repr__(self):
         repr_str = f"{colors.BLUE}{self._repr_type()}(len={colors.MINT}{formats.UNDERLINE}{len(self)}{formats.RESET}{colors.BLUE}, ["
 
diff --git a/genesis/engine/entities/mpm_entity.py b/genesis/engine/entities/mpm_entity.py
index 104b4957..d9a8a4ce 100644
--- a/genesis/engine/entities/mpm_entity.py
+++ b/genesis/engine/entities/mpm_entity.py
@@ -185,6 +185,27 @@ def set_muscle(self, muscle_group=None, muscle_direction=None):
 
             self.set_muscle_direction(muscle_direction)
 
+    def set_free(self, free):
+        self._assert_active()
+
+        self.solver._kernel_set_free(
+            self._particle_start,
+            self._n_particles,
+            free,
+        )
+
+    def get_free(self):
+        self._assert_active()
+
+        free = gs.zeros((self._n_particles,), dtype=int, requires_grad=False, scene=self._scene)
+        self.solver._kernel_get_free(
+            self._particle_start,
+            self._n_particles,
+            free,
+        )
+
+        return free
+
     @ti.kernel
     def clear_grad(self, f: ti.i32):
         for i in range(self.n_particles):
diff --git a/genesis/engine/entities/rigid_entity/rigid_entity.py b/genesis/engine/entities/rigid_entity/rigid_entity.py
index 0517c80a..955af57c 100644
--- a/genesis/engine/entities/rigid_entity/rigid_entity.py
+++ b/genesis/engine/entities/rigid_entity/rigid_entity.py
@@ -666,7 +666,8 @@ def _func_get_jacobian(self, tgt_link_idx, i_b, pos_mask, rot_mask):
         tgt_link_pos = self._solver.links_state[tgt_link_idx, i_b].pos
         i_l = tgt_link_idx
         while i_l > -1:
-            l_info = self._solver.links_info[i_l]
+            I_l = [i_l, i_b] if ti.static(self.solver._options.batch_links_info) else i_l
+            l_info = self._solver.links_info[I_l]
             l_state = self._solver.links_state[i_l, i_b]
 
             if l_info.joint_type == gs.JOINT_TYPE.FIXED:
@@ -674,8 +675,9 @@ def _func_get_jacobian(self, tgt_link_idx, i_b, pos_mask, rot_mask):
 
             elif l_info.joint_type == gs.JOINT_TYPE.REVOLUTE:
                 i_d = l_info.dof_start
+                I_d = [i_d, i_b] if ti.static(self.solver._options.batch_dofs_info) else i_d
                 i_d_jac = i_d - self._dof_start
-                rotation = gu.ti_transform_by_quat(self._solver.dofs_info[i_d].motion_ang, l_state.quat)
+                rotation = gu.ti_transform_by_quat(self._solver.dofs_info[I_d].motion_ang, l_state.quat)
                 translation = rotation.cross(tgt_link_pos - l_state.pos)
 
                 self._jacobian[0, i_d_jac, i_b] = translation[0] * pos_mask[0]
@@ -687,8 +689,9 @@ def _func_get_jacobian(self, tgt_link_idx, i_b, pos_mask, rot_mask):
 
             elif l_info.joint_type == gs.JOINT_TYPE.PRISMATIC:
                 i_d = l_info.dof_start
+                I_d = [i_d, i_b] if ti.static(self.solver._options.batch_dofs_info) else i_d
                 i_d_jac = i_d - self._dof_start
-                translation = gu.ti_transform_by_quat(self._solver.dofs_info[i_d].motion_vel, l_state.quat)
+                translation = gu.ti_transform_by_quat(self._solver.dofs_info[I_d].motion_vel, l_state.quat)
 
                 self._jacobian[0, i_d_jac, i_b] = translation[0] * pos_mask[0]
                 self._jacobian[1, i_d_jac, i_b] = translation[1] * pos_mask[1]
@@ -706,7 +709,8 @@ def _func_get_jacobian(self, tgt_link_idx, i_b, pos_mask, rot_mask):
                 for i_d_ in range(3):
                     i_d = l_info.dof_start + i_d_ + 3
                     i_d_jac = i_d - self._dof_start
-                    rotation = self._solver.dofs_info[i_d].motion_ang
+                    I_d = [i_d, i_b] if ti.static(self.solver._options.batch_dofs_info) else i_d
+                    rotation = self._solver.dofs_info[I_d].motion_ang
                     translation = rotation.cross(tgt_link_pos - l_state.pos)
 
                     self._jacobian[0, i_d_jac, i_b] = translation[0] * pos_mask[0]
@@ -1156,8 +1160,14 @@ def _kernel_inverse_kinematics(
                     # Resample init q
                     if respect_joint_limit and i_sample < max_samples - 1:
                         for i_l in range(self.link_start, self.link_end):
-                            l_info = self._solver.links_info[i_l]
-                            dof_info = self._solver.dofs_info[l_info.dof_start]
+                            I_l = [i_l, i_b] if ti.static(self.solver._options.batch_links_info) else i_l
+                            l_info = self._solver.links_info[I_l]
+                            I_dof_start = (
+                                [l_info.dof_start, i_b]
+                                if ti.static(self.solver._options.batch_dofs_info)
+                                else l_info.dof_start
+                            )
+                            dof_info = self._solver.dofs_info[I_dof_start]
                             q_start = l_info.q_start
 
                             if l_info.joint_type == gs.JOINT_TYPE.FREE:
@@ -1561,6 +1571,14 @@ def get_links_ang(self, envs_idx=None):
         """
         return self._solver.get_links_ang(np.arange(self.link_start, self.link_end), envs_idx)
 
+    @gs.assert_built
+    def get_links_inertial_mass(self, ls_idx_local=None, envs_idx=None):
+        return self._solver.get_links_inertial_mass(self._get_ls_idx(ls_idx_local), envs_idx)
+
+    @gs.assert_built
+    def get_links_invweight(self, ls_idx_local=None, envs_idx=None):
+        return self._solver.get_links_invweight(self._get_ls_idx(ls_idx_local), envs_idx)
+
     @gs.assert_built
     def set_pos(self, pos, zero_velocity=True, envs_idx=None):
         """
@@ -1694,6 +1712,18 @@ def _get_dofs_idx_local(self, dofs_idx_local=None):
     def _get_dofs_idx(self, dofs_idx_local=None):
         return self._get_dofs_idx_local(dofs_idx_local) + self._dof_start
 
+    def _get_ls_idx_local(self, ls_idx_local=None):
+        if ls_idx_local is None:
+            ls_idx_local = torch.arange(self.n_links, dtype=torch.int32, device=gs.device)
+        else:
+            ls_idx_local = torch.as_tensor(ls_idx_local, dtype=gs.tc_int)
+            if (ls_idx_local < 0).any() or (ls_idx_local >= self.n_links).any():
+                gs.raise_exception("`ls_idx_local` exceeds valid range.")
+        return ls_idx_local
+
+    def _get_ls_idx(self, ls_idx_local=None):
+        return self._get_ls_idx_local(ls_idx_local) + self._link_start
+
     @gs.assert_built
     def set_qpos(self, qpos, qs_idx_local=None, zero_velocity=True, envs_idx=None):
         """
@@ -1716,7 +1746,7 @@ def set_qpos(self, qpos, qs_idx_local=None, zero_velocity=True, envs_idx=None):
             self.zero_all_dofs_velocity(envs_idx)
 
     @gs.assert_built
-    def set_dofs_kp(self, kp, dofs_idx_local=None):
+    def set_dofs_kp(self, kp, dofs_idx_local=None, envs_idx=None):
         """
         Set the entity's dofs' positional gains for the PD controller.
 
@@ -1726,12 +1756,14 @@ def set_dofs_kp(self, kp, dofs_idx_local=None):
             The positional gains to set.
         dofs_idx_local : None | array_like, optional
             The indices of the dofs to set. If None, all dofs will be set. Note that here this uses the local `q_idx`, not the scene-level one. Defaults to None.
+        envs_idx : None | array_like, optional
+            The indices of the environments. If None, all environments will be considered. Defaults to None.
         """
 
-        self._solver.set_dofs_kp(kp, self._get_dofs_idx(dofs_idx_local))
+        self._solver.set_dofs_kp(kp, self._get_dofs_idx(dofs_idx_local), envs_idx)
 
     @gs.assert_built
-    def set_dofs_kv(self, kv, dofs_idx_local=None):
+    def set_dofs_kv(self, kv, dofs_idx_local=None, envs_idx=None):
         """
         Set the entity's dofs' velocity gains for the PD controller.
 
@@ -1741,11 +1773,13 @@ def set_dofs_kv(self, kv, dofs_idx_local=None):
             The velocity gains to set.
         dofs_idx_local : None | array_like, optional
             The indices of the dofs to set. If None, all dofs will be set. Note that here this uses the local `q_idx`, not the scene-level one. Defaults to None.
+        envs_idx : None | array_like, optional
+            The indices of the environments. If None, all environments will be considered. Defaults to None.
         """
-        self._solver.set_dofs_kv(kv, self._get_dofs_idx(dofs_idx_local))
+        self._solver.set_dofs_kv(kv, self._get_dofs_idx(dofs_idx_local), envs_idx)
 
     @gs.assert_built
-    def set_dofs_force_range(self, lower, upper, dofs_idx_local=None):
+    def set_dofs_force_range(self, lower, upper, dofs_idx_local=None, envs_idx=None):
         """
         Set the entity's dofs' force range.
 
@@ -1757,9 +1791,27 @@ def set_dofs_force_range(self, lower, upper, dofs_idx_local=None):
             The upper bounds of the force range.
         dofs_idx_local : None | array_like, optional
             The indices of the dofs to set. If None, all dofs will be set. Note that here this uses the local `q_idx`, not the scene-level one. Defaults to None.
+        envs_idx : None | array_like, optional
+            The indices of the environments. If None, all environments will be considered. Defaults to None.
         """
 
-        self._solver.set_dofs_force_range(lower, upper, self._get_dofs_idx(dofs_idx_local))
+        self._solver.set_dofs_force_range(lower, upper, self._get_dofs_idx(dofs_idx_local), envs_idx)
+
+    @gs.assert_built
+    def set_dofs_stiffness(self, stiffness, dofs_idx_local=None, envs_idx=None):
+        self._solver.set_dofs_stiffness(stiffness, self._get_dofs_idx(dofs_idx_local), envs_idx)
+
+    @gs.assert_built
+    def set_dofs_invweight(self, invweight, dofs_idx_local=None, envs_idx=None):
+        self._solver.set_dofs_invweight(invweight, self._get_dofs_idx(dofs_idx_local), envs_idx)
+
+    @gs.assert_built
+    def set_dofs_armature(self, armature, dofs_idx_local=None, envs_idx=None):
+        self._solver.set_dofs_armature(armature, self._get_dofs_idx(dofs_idx_local), envs_idx)
+
+    @gs.assert_built
+    def set_dofs_damping(self, damping, dofs_idx_local=None, envs_idx=None):
+        self._solver.set_dofs_damping(damping, self._get_dofs_idx(dofs_idx_local), envs_idx)
 
     @gs.assert_built
     def set_dofs_velocity(self, velocity, dofs_idx_local=None, envs_idx=None):
@@ -1948,7 +2000,7 @@ def get_dofs_position(self, dofs_idx_local=None, envs_idx=None):
         return self._solver.get_dofs_position(self._get_dofs_idx(dofs_idx_local), envs_idx)
 
     @gs.assert_built
-    def get_dofs_kp(self, dofs_idx_local=None):
+    def get_dofs_kp(self, dofs_idx_local=None, envs_idx=None):
         """
         Get the positional gain (kp) for the entity's dofs used by the PD controller.
 
@@ -1956,16 +2008,18 @@ def get_dofs_kp(self, dofs_idx_local=None):
         ----------
         dofs_idx_local : None | array_like, optional
             The indices of the dofs to get. If None, all dofs will be returned. Note that here this uses the local `q_idx`, not the scene-level one. Defaults to None.
+        envs_idx : None | array_like, optional
+            The indices of the environments. If None, all environments will be considered. Defaults to None.
 
         Returns
         -------
-        kp : torch.Tensor, shape (n_dofs,)
+        kp : torch.Tensor, shape (n_dofs,) or (n_envs, n_dofs)
             The positional gain (kp) for the entity's dofs.
         """
-        return self._solver.get_dofs_kp(self._get_dofs_idx(dofs_idx_local))
+        return self._solver.get_dofs_kp(self._get_dofs_idx(dofs_idx_local), envs_idx)
 
     @gs.assert_built
-    def get_dofs_kv(self, dofs_idx_local=None):
+    def get_dofs_kv(self, dofs_idx_local=None, envs_idx=None):
         """
         Get the velocity gain (kv) for the entity's dofs used by the PD controller.
 
@@ -1973,16 +2027,18 @@ def get_dofs_kv(self, dofs_idx_local=None):
         ----------
         dofs_idx_local : None | array_like, optional
             The indices of the dofs to get. If None, all dofs will be returned. Note that here this uses the local `q_idx`, not the scene-level one. Defaults to None.
+        envs_idx : None | array_like, optional
+            The indices of the environments. If None, all environments will be considered. Defaults to None.
 
         Returns
         -------
-        kv : torch.Tensor, shape (n_dofs,)
+        kv : torch.Tensor, shape (n_dofs,) or (n_envs, n_dofs)
             The velocity gain (kv) for the entity's dofs.
         """
-        return self._solver.get_dofs_kv(self._get_dofs_idx(dofs_idx_local))
+        return self._solver.get_dofs_kv(self._get_dofs_idx(dofs_idx_local), envs_idx)
 
     @gs.assert_built
-    def get_dofs_force_range(self, dofs_idx_local=None):
+    def get_dofs_force_range(self, dofs_idx_local=None, envs_idx=None):
         """
         Get the force range (min and max limits) for the entity's dofs.
 
@@ -1990,18 +2046,20 @@ def get_dofs_force_range(self, dofs_idx_local=None):
         ----------
         dofs_idx_local : None | array_like, optional
             The indices of the dofs to get. If None, all dofs will be returned. Note that here this uses the local `q_idx`, not the scene-level one. Defaults to None.
+        envs_idx : None | array_like, optional
+            The indices of the environments. If None, all environments will be considered. Defaults to None.
 
         Returns
         -------
-        lower_limit : torch.Tensor, shape (n_dofs,)
+        lower_limit : torch.Tensor, shape (n_dofs,) or (n_envs, n_dofs)
             The lower limit of the force range for the entity's dofs.
-        upper_limit : torch.Tensor, shape (n_dofs,)
+        upper_limit : torch.Tensor, shape (n_dofs,) or (n_envs, n_dofs)
             The upper limit of the force range for the entity's dofs.
         """
-        return self._solver.get_dofs_force_range(self._get_dofs_idx(dofs_idx_local))
+        return self._solver.get_dofs_force_range(self._get_dofs_idx(dofs_idx_local), envs_idx)
 
     @gs.assert_built
-    def get_dofs_limit(self, dofs_idx=None):
+    def get_dofs_limit(self, dofs_idx=None, envs_idx=None):
         """
         Get the positional limits (min and max) for the entity's dofs.
 
@@ -2009,15 +2067,33 @@ def get_dofs_limit(self, dofs_idx=None):
         ----------
         dofs_idx : None | array_like, optional
             The indices of the dofs to get. If None, all dofs will be returned. Note that here this uses the local `q_idx`, not the scene-level one. Defaults to None.
+        envs_idx : None | array_like, optional
+            The indices of the environments. If None, all environments will be considered. Defaults to None.
 
         Returns
         -------
-        lower_limit : torch.Tensor, shape (n_dofs,)
+        lower_limit : torch.Tensor, shape (n_dofs,) or (n_envs, n_dofs)
             The lower limit of the positional limits for the entity's dofs.
-        upper_limit : torch.Tensor, shape (n_dofs,)
+        upper_limit : torch.Tensor, shape (n_dofs,) or (n_envs, n_dofs)
             The upper limit of the positional limits for the entity's dofs.
         """
-        return self._solver.get_dofs_limit(self._get_dofs_idx(dofs_idx))
+        return self._solver.get_dofs_limit(self._get_dofs_idx(dofs_idx), envs_idx)
+
+    @gs.assert_built
+    def get_dofs_stiffness(self, dofs_idx_local=None, envs_idx=None):
+        return self._solver.get_dofs_stiffness(self._get_dofs_idx(dofs_idx_local), envs_idx)
+
+    @gs.assert_built
+    def get_dofs_invweight(self, dofs_idx_local=None, envs_idx=None):
+        return self._solver.get_dofs_invweight(self._get_dofs_idx(dofs_idx_local), envs_idx)
+
+    @gs.assert_built
+    def get_dofs_armature(self, dofs_idx_local=None, envs_idx=None):
+        return self._solver.get_dofs_armature(self._get_dofs_idx(dofs_idx_local), envs_idx)
+
+    @gs.assert_built
+    def get_dofs_damping(self, dofs_idx_local=None, envs_idx=None):
+        return self._solver.get_dofs_damping(self._get_dofs_idx(dofs_idx_local), envs_idx)
 
     @gs.assert_built
     def zero_all_dofs_velocity(self, envs_idx=None):
@@ -2251,6 +2327,14 @@ def set_COM_shift(self, com_shift, link_indices, envs_idx=None):
             link_indices[i] += self._link_start
         self._solver.set_links_COM_shift(com_shift, link_indices, envs_idx)
 
+    @gs.assert_built
+    def set_links_inertial_mass(self, inertial_mass, ls_idx_local=None, envs_idx=None):
+        self._solver.set_links_inertial_mass(inertial_mass, self._get_ls_idx(ls_idx_local), envs_idx)
+
+    @gs.assert_built
+    def set_links_invweight(self, invweight, ls_idx_local=None, envs_idx=None):
+        self._solver.set_links_invweight(invweight, self._get_ls_idx(ls_idx_local), envs_idx)
+
     @gs.assert_built
     def get_mass(self):
         """
diff --git a/genesis/engine/entities/rigid_entity/rigid_joint.py b/genesis/engine/entities/rigid_entity/rigid_joint.py
index 660971d9..e81c25b2 100644
--- a/genesis/engine/entities/rigid_entity/rigid_joint.py
+++ b/genesis/engine/entities/rigid_entity/rigid_joint.py
@@ -90,7 +90,8 @@ def get_pos(self):
     def _kernel_get_pos(self, tensor: ti.types.ndarray()):
 
         for i_b in range(self._solver._B):
-            l_info = self._solver.links_info[self._idx, i_b]
+            I_l = [self._idx, i_b] if ti.static(self._solver._options.batch_links_info) else self._idx
+            l_info = self._solver.links_info[I_l]
             i_p = l_info.parent_idx
 
             p_pos = ti.Vector.zero(gs.ti_float, 3)
@@ -124,7 +125,8 @@ def get_quat(self):
     def _kernel_get_quat(self, tensor: ti.types.ndarray()):
 
         for i_b in range(self._solver._B):
-            l_info = self._solver.links_info[self._idx, i_b]
+            I_l = [self._idx, i_b] if ti.static(self._solver._options.batch_links_info) else self._idx
+            l_info = self._solver.links_info[I_l]
             i_p = l_info.parent_idx
 
             p_pos = ti.Vector.zero(gs.ti_float, 3)
diff --git a/genesis/engine/solvers/mpm_solver.py b/genesis/engine/solvers/mpm_solver.py
index 6356a76b..8d490022 100644
--- a/genesis/engine/solvers/mpm_solver.py
+++ b/genesis/engine/solvers/mpm_solver.py
@@ -93,6 +93,7 @@ def init_particle_fields(self):
             mat_idx=gs.ti_int,
             mass=gs.ti_float,
             default_Jp=gs.ti_float,
+            free=gs.ti_int,
             # for muscle
             muscle_group=gs.ti_int,
             muscle_direction=gs.ti_vec3,
@@ -384,6 +385,9 @@ def p2g(self, f: ti.i32):
                         )
                         self.grid[f, base - self._grid_offset + offset].mass += weight * self.particles_info[i].mass
 
+                    if self.particles_info[i].free == 0:  # non-free particles behave as boundary conditions
+                        self.grid[f, base - self._grid_offset + offset].vel_in = ti.Vector.zero(gs.ti_float, 3)
+
     @ti.kernel
     def g2p(self, f: ti.i32):
         for i in range(self._n_particles):
@@ -658,6 +662,7 @@ def _kernel_add_particles(
             self.particles_info[i_global].mat_idx = mat_idx
             self.particles_info[i_global].default_Jp = mat_default_Jp
             self.particles_info[i_global].mass = self._p_vol * mat_rho
+            self.particles_info[i_global].free = 1
             self.particles_info[i_global].muscle_group = 0
             self.particles_info[i_global].muscle_direction = ti.Vector([0.0, 0.0, 1.0], dt=gs.ti_float)
 
@@ -816,6 +821,28 @@ def _kernel_set_muscle_direction(
             for j in ti.static(range(3)):
                 self.particles_info[i_global].muscle_direction[j] = muscle_direction[i, j]
 
+    @ti.kernel
+    def _kernel_set_free(
+        self,
+        particle_start: ti.i32,
+        n_particles: ti.i32,
+        free: ti.types.ndarray(),
+    ):
+        for i in range(n_particles):
+            i_global = i + particle_start
+            self.particles_info[i_global].free = free[i]
+
+    @ti.kernel
+    def _kernel_get_free(
+        self,
+        particle_start: ti.i32,
+        n_particles: ti.i32,
+        free: ti.types.ndarray(),
+    ):
+        for i in range(n_particles):
+            i_global = i + particle_start
+            free[i] = self.particles_info[i_global].free
+
     @ti.kernel
     def _kernel_get_state(
         self,
diff --git a/genesis/engine/solvers/rigid/collider_decomp.py b/genesis/engine/solvers/rigid/collider_decomp.py
index 968ad214..35a7756f 100644
--- a/genesis/engine/solvers/rigid/collider_decomp.py
+++ b/genesis/engine/solvers/rigid/collider_decomp.py
@@ -51,6 +51,10 @@ def _init_collision_fields(self):
         links_root_idx = self._solver.links_info.root_idx.to_numpy()
         links_parent_idx = self._solver.links_info.parent_idx.to_numpy()
         links_is_fixed = self._solver.links_info.is_fixed.to_numpy()
+        if self._solver._options.batch_links_info:
+            links_root_idx = links_root_idx[:, 0]
+            links_parent_idx = links_parent_idx[:, 0]
+            links_is_fixed = links_is_fixed[:, 0]
         n_possible_pairs = 0
         for i in range(self._solver.n_geoms):
             for j in range(i + 1, self._solver.n_geoms):
@@ -164,10 +168,13 @@ def clear(self):
                     i_la = self.contact_data[i_c, i_b].link_a
                     i_lb = self.contact_data[i_c, i_b].link_b
 
+                    I_la = [i_la, i_b] if ti.static(self._solver._options.batch_links_info) else i_la
+                    I_lb = [i_lb, i_b] if ti.static(self._solver._options.batch_links_info) else i_lb
+
                     # pair of hibernated-fixed links -> hibernated contact
                     # TODO: we should also include hibernated-hibernated links and wake up the whole contact island once a new collision is detected
-                    if (self._solver.links_state[i_la, i_b].hibernated and self._solver.links_info[i_lb].is_fixed) or (
-                        self._solver.links_state[i_lb, i_b].hibernated and self._solver.links_info[i_la].is_fixed
+                    if (self._solver.links_state[i_la, i_b].hibernated and self._solver.links_info[I_lb].is_fixed) or (
+                        self._solver.links_state[i_lb, i_b].hibernated and self._solver.links_info[I_la].is_fixed
                     ):
                         i_c_hibernated = self.n_contacts_hibernated[i_b]
                         if i_c != i_c_hibernated:
@@ -549,6 +556,8 @@ def _func_update_aabbs(self):
     def _func_check_collision_valid(self, i_ga, i_gb, i_b):
         i_la = self._solver.geoms_info[i_ga].link_idx
         i_lb = self._solver.geoms_info[i_gb].link_idx
+        I_la = [i_la, i_b] if ti.static(self._solver._options.batch_links_info) else i_la
+        I_lb = [i_lb, i_b] if ti.static(self._solver._options.batch_links_info) else i_lb
         is_valid = True
 
         # geoms in the same link
@@ -558,22 +567,22 @@ def _func_check_collision_valid(self, i_ga, i_gb, i_b):
         # self collision
         if (
             ti.static(not self._solver._enable_self_collision)
-            and self._solver.links_info[i_la].root_idx == self._solver.links_info[i_lb].root_idx
+            and self._solver.links_info[I_la].root_idx == self._solver.links_info[I_lb].root_idx
         ):
             is_valid = False
 
         # adjacent links
-        if self._solver.links_info[i_la].parent_idx == i_lb or self._solver.links_info[i_lb].parent_idx == i_la:
+        if self._solver.links_info[I_la].parent_idx == i_lb or self._solver.links_info[I_lb].parent_idx == i_la:
             is_valid = False
 
         # pair of fixed links
-        if self._solver.links_info[i_la].is_fixed and self._solver.links_info[i_lb].is_fixed:
+        if self._solver.links_info[I_la].is_fixed and self._solver.links_info[I_lb].is_fixed:
             is_valid = False
 
         # hibernated <-> fixed links
         if ti.static(self._solver._use_hibernation):
-            if (self._solver.links_state[i_la, i_b].hibernated and self._solver.links_info[i_lb].is_fixed) or (
-                self._solver.links_state[i_lb, i_b].hibernated and self._solver.links_info[i_la].is_fixed
+            if (self._solver.links_state[i_la, i_b].hibernated and self._solver.links_info[I_lb].is_fixed) or (
+                self._solver.links_state[i_lb, i_b].hibernated and self._solver.links_info[I_la].is_fixed
             ):
                 is_valid = False
 
@@ -994,7 +1003,9 @@ def _func_mpr(self, i_ga, i_gb, i_b):
 
         i_la = self._solver.geoms_info[i_ga].link_idx
         i_lb = self._solver.geoms_info[i_gb].link_idx
-        is_self_pair = self._solver.links_info.root_idx[i_la] == self._solver.links_info.root_idx[i_lb]
+        I_la = [i_la, i_b] if ti.static(self._solver._options.batch_links_info) else i_la
+        I_lb = [i_lb, i_b] if ti.static(self._solver._options.batch_links_info) else i_lb
+        is_self_pair = self._solver.links_info.root_idx[I_la] == self._solver.links_info.root_idx[I_lb]
         multi_contact = (
             self._solver.geoms_info[i_ga].type != gs.GEOM_TYPE.SPHERE
             and self._solver.geoms_info[i_gb].type != gs.GEOM_TYPE.SPHERE
diff --git a/genesis/engine/solvers/rigid/constraint_solver_decomp.py b/genesis/engine/solvers/rigid/constraint_solver_decomp.py
index c680c31f..d3592b92 100644
--- a/genesis/engine/solvers/rigid/constraint_solver_decomp.py
+++ b/genesis/engine/solvers/rigid/constraint_solver_decomp.py
@@ -105,14 +105,16 @@ def add_collision_constraints(self):
                 impact = self._collider.contact_data[i_col, i_b]
                 link_a = impact.link_a
                 link_b = impact.link_b
+                link_a_maybe_batch = [link_a, i_b] if ti.static(self._solver._options.batch_links_info) else link_a
+                link_b_maybe_batch = [link_b, i_b] if ti.static(self._solver._options.batch_links_info) else link_b
                 f = impact.friction
                 pos = impact.pos
 
                 d1, d2 = gu.orthogonals(impact.normal)
 
-                t = self._solver.links_info[link_a].invweight + self._solver.links_info[link_b].invweight * (
-                    link_b > -1
-                )
+                t = self._solver.links_info[link_a_maybe_batch].invweight + self._solver.links_info[
+                    link_b_maybe_batch
+                ].invweight * (link_b > -1)
                 for i in range(4):
                     n = -d1 * f - impact.normal
                     if i == 1:
@@ -143,10 +145,13 @@ def add_collision_constraints(self):
                             link = link_b
 
                         while link > -1:
+                            link_maybe_batch = (
+                                [link, i_b] if ti.static(self._solver._options.batch_links_info) else link
+                            )
 
                             # reverse order to make sure dofs in each row of self.jac_relevant_dofs is strictly descending
-                            for i_d_ in range(self._solver.links_info[link].n_dofs):
-                                i_d = self._solver.links_info[link].dof_end - 1 - i_d_
+                            for i_d_ in range(self._solver.links_info[link_maybe_batch].n_dofs):
+                                i_d = self._solver.links_info[link_maybe_batch].dof_end - 1 - i_d_
 
                                 cdof_ang = self._solver.dofs_state[i_d, i_b].cdof_ang
                                 cdot_vel = self._solver.dofs_state[i_d, i_b].cdof_vel
@@ -164,7 +169,7 @@ def add_collision_constraints(self):
                                     self.jac_relevant_dofs[n_con, con_n_relevant_dofs, i_b] = i_d
                                     con_n_relevant_dofs += 1
 
-                            link = self._solver.links_info[link].parent_idx
+                            link = self._solver.links_info[link_maybe_batch].parent_idx
 
                     if ti.static(self.sparse_solve):
                         self.jac_n_relevant_dofs[n_con, i_b] = con_n_relevant_dofs
@@ -183,21 +188,23 @@ def add_joint_limit_constraints(self):
         ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.PARTIAL)
         for i_b in range(self._B):
             for i_l in range(self._solver.n_links):
-                l_info = self._solver.links_info[i_l]
+                I_l = [i_l, i_b] if ti.static(self._solver._options.batch_links_info) else i_l
+                l_info = self._solver.links_info[I_l]
                 if l_info.joint_type == gs.JOINT_TYPE.REVOLUTE or l_info.joint_type == gs.JOINT_TYPE.PRISMATIC:
 
                     i_q = l_info.q_start
                     i_d = l_info.dof_start
-                    pos_min = self._solver.qpos[i_q, i_b] - self._solver.dofs_info[i_d].limit[0]
-                    pos_max = self._solver.dofs_info[i_d].limit[1] - self._solver.qpos[i_q, i_b]
+                    I_d = [i_d, i_b] if ti.static(self._solver._options.batch_dofs_info) else i_d
+                    pos_min = self._solver.qpos[i_q, i_b] - self._solver.dofs_info[I_d].limit[0]
+                    pos_max = self._solver.dofs_info[I_d].limit[1] - self._solver.qpos[i_q, i_b]
                     pos = min(min(pos_min, pos_max), 0)
 
                     side = ((pos_min < pos_max) * 2 - 1) * (pos < 0)
 
                     jac = side
                     jac_qvel = jac * self._solver.dofs_state[i_d, i_b].vel
-                    imp, aref = gu.imp_aref(self._solver.dofs_info[i_d].sol_params, pos, jac_qvel)
-                    diag = self._solver.dofs_info[i_d].invweight * (pos < 0) * (1 - imp) / (imp + gs.EPS)
+                    imp, aref = gu.imp_aref(self._solver.dofs_info[I_d].sol_params, pos, jac_qvel)
+                    diag = self._solver.dofs_info[I_d].invweight * (pos < 0) * (1 - imp) / (imp + gs.EPS)
                     aref = aref * (pos < 0)
                     if pos < 0:
                         n_con = self.n_constraints[i_b]
diff --git a/genesis/engine/solvers/rigid/constraint_solver_decomp_island.py b/genesis/engine/solvers/rigid/constraint_solver_decomp_island.py
index e33867a0..4a7d2257 100644
--- a/genesis/engine/solvers/rigid/constraint_solver_decomp_island.py
+++ b/genesis/engine/solvers/rigid/constraint_solver_decomp_island.py
@@ -131,12 +131,16 @@ def add_collision_constraints(self, island, i_b):
             impact = self._collider.contact_data[i_col, i_b]
             link_a = impact.link_a
             link_b = impact.link_b
+            link_a_maybe_batch = [link_a, i_b] if ti.static(self._solver._options.batch_links_info) else link_a
+            link_b_maybe_batch = [link_b, i_b] if ti.static(self._solver._options.batch_links_info) else link_b
             f = impact.friction
             pos = impact.pos
 
             d1, d2 = gu.orthogonals(impact.normal)
 
-            t = self._solver.links_info[link_a].invweight + self._solver.links_info[link_b].invweight * (link_b > -1)
+            t = self._solver.links_info[link_a_maybe_batch].invweight + self._solver.links_info[
+                link_b_maybe_batch
+            ].invweight * (link_b > -1)
 
             for i in range(4):
                 n = -d1 * f - impact.normal
@@ -170,10 +174,11 @@ def add_collision_constraints(self, island, i_b):
                         link = link_b
 
                     while link > -1:
+                        link_maybe_batch = [link, i_b] if ti.static(self._solver._options.batch_links_info) else link
 
                         # reverse order to make sure dofs in each row of self.jac_relevant_dofs is strictly descending
                         for i_d_ in range(self._solver.links_info[link].n_dofs):
-                            i_d = self._solver.links_info[link].dof_end - 1 - i_d_
+                            i_d = self._solver.links_info[link_maybe_batch].dof_end - 1 - i_d_
 
                             cdof_ang = self._solver.dofs_state[i_d, i_b].cdof_ang
                             cdot_vel = self._solver.dofs_state[i_d, i_b].cdof_vel
@@ -190,7 +195,7 @@ def add_collision_constraints(self, island, i_b):
                                 self.jac_relevant_dofs[n_con, con_n_relevant_dofs, i_b] = i_d
                                 con_n_relevant_dofs += 1
 
-                        link = self._solver.links_info[link].parent_idx
+                        link = self._solver.links_info[link_maybe_batch].parent_idx
 
                 if ti.static(self.sparse_solve):
                     self.jac_n_relevant_dofs[n_con, i_b] = con_n_relevant_dofs
@@ -207,8 +212,8 @@ def add_collision_constraints(self, island, i_b):
 
             if ti.static(self._solver._use_hibernation):
                 # wake up entities
-                self._solver._func_wakeup_entity(self._solver.links_info[link_a].entity_idx, i_b)
-                self._solver._func_wakeup_entity(self._solver.links_info[link_b].entity_idx, i_b)
+                self._solver._func_wakeup_entity(self._solver.links_info[link_a_maybe_batch].entity_idx, i_b)
+                self._solver._func_wakeup_entity(self._solver.links_info[link_b_maybe_batch].entity_idx, i_b)
 
     @ti.func
     def add_joint_limit_constraints(self, island, i_b):
@@ -219,22 +224,23 @@ def add_joint_limit_constraints(self, island, i_b):
             e_info = self.entities_info[i_e]
 
             for i_l in range(e_info.link_start, e_info.link_end):
-
-                l_info = self._solver.links_info[i_l]
+                I_l = [i_l, i_b] if ti.static(self._solver._options.batch_links_info) else i_l
+                l_info = self._solver.links_info[I_l]
                 if l_info.joint_type == gs.JOINT_TYPE.REVOLUTE or l_info.joint_type == gs.JOINT_TYPE.PRISMATIC:
 
                     i_q = l_info.q_start
                     i_d = l_info.dof_start
-                    pos_min = self._solver.qpos[i_q, i_b] - self._solver.dofs_info[i_d].limit[0]
-                    pos_max = self._solver.dofs_info[i_d].limit[1] - self._solver.qpos[i_q, i_b]
+                    I_d = [i_d, i_b] if ti.static(self._solver._options.batch_dofs_info) else i_d
+                    pos_min = self._solver.qpos[i_q, i_b] - self._solver.dofs_info[I_d].limit[0]
+                    pos_max = self._solver.dofs_info[I_d].limit[1] - self._solver.qpos[i_q, i_b]
                     pos = min(min(pos_min, pos_max), 0)
 
                     side = ((pos_min < pos_max) * 2 - 1) * (pos < 0)
 
                     jac = side
                     jac_qvel = jac * self._solver.dofs_state[i_d, i_b].vel
-                    imp, aref = gu.imp_aref(self._solver.dofs_info[i_d].sol_params, pos, jac_qvel)
-                    diag = self._solver.dofs_info[i_d].invweight * (pos < 0) * (1 - imp) / (imp + gs.EPS)
+                    imp, aref = gu.imp_aref(self._solver.dofs_info[I_d].sol_params, pos, jac_qvel)
+                    diag = self._solver.dofs_info[I_d].invweight * (pos < 0) * (1 - imp) / (imp + gs.EPS)
                     aref = aref * (pos < 0)
                     if pos < 0:
                         n_con = self.n_constraints[i_b]
diff --git a/genesis/engine/solvers/rigid/contact_island.py b/genesis/engine/solvers/rigid/contact_island.py
index 2914b503..97ae3b93 100644
--- a/genesis/engine/solvers/rigid/contact_island.py
+++ b/genesis/engine/solvers/rigid/contact_island.py
@@ -66,8 +66,11 @@ def clear(self):
 
     @ti.func
     def add_edge(self, link_a, link_b, i_b):
-        ea = self.solver.links_info[link_a].entity_idx
-        eb = self.solver.links_info[link_b].entity_idx
+        link_a_maybe_batch = [link_a, i_b] if ti.static(self.solver._options.batch_links_info) else link_a
+        link_b_maybe_batch = [link_b, i_b] if ti.static(self.solver._options.batch_links_info) else link_b
+
+        ea = self.solver.links_info[link_a_maybe_batch].entity_idx
+        eb = self.solver.links_info[link_b_maybe_batch].entity_idx
 
         self.entity_edge[ea, i_b].n = self.entity_edge[ea, i_b].n + 1
         self.entity_edge[eb, i_b].n = self.entity_edge[eb, i_b].n + 1
@@ -100,9 +103,11 @@ def postprocess_island(self):
                 impact = self.collider.contact_data[i_col, i_b]
                 link_a = impact.link_a
                 link_b = impact.link_b
+                link_a_maybe_batch = [link_a, i_b] if ti.static(self.solver._options.batch_links_info) else link_a
+                link_b_maybe_batch = [link_b, i_b] if ti.static(self.solver._options.batch_links_info) else link_b
 
-                ea = self.solver.links_info[link_a].entity_idx
-                eb = self.solver.links_info[link_b].entity_idx
+                ea = self.solver.links_info[link_a_maybe_batch].entity_idx
+                eb = self.solver.links_info[link_b_maybe_batch].entity_idx
 
                 island_a = self.entity_island[ea, i_b]
                 island_b = self.entity_island[eb, i_b]
diff --git a/genesis/engine/solvers/rigid/rigid_solver_decomp.py b/genesis/engine/solvers/rigid/rigid_solver_decomp.py
index a6abeb53..b6dd3a17 100644
--- a/genesis/engine/solvers/rigid/rigid_solver_decomp.py
+++ b/genesis/engine/solvers/rigid/rigid_solver_decomp.py
@@ -159,6 +159,11 @@ def _init_invweight(self):
         dof_end = self.links_info.dof_end.to_numpy()
         n_dofs = self.links_info.n_dofs.to_numpy()
         parent_idx = self.links_info.parent_idx.to_numpy()
+        if self._options.batch_links_info:
+            dof_start = dof_start[:, 0]
+            dof_end = dof_end[:, 0]
+            n_dofs = n_dofs[:, 0]
+            parent_idx = parent_idx[:, 0]
 
         offsets = self.links_state.i_pos.to_numpy()[:, 0, :]
 
@@ -209,9 +214,9 @@ def _kernel_init_invweight(
         invweight: ti.types.ndarray(),
     ):
         ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.PARTIAL)
-        for i, i_b in ti.ndrange(self.n_links, self._B):
-            if self.links_info[i].invweight < 0:
-                self.links_info[i].invweight = invweight[i]
+        for I in ti.grouped(self.links_info):
+            if self.links_info[I].invweight < 0:
+                self.links_info[I].invweight = invweight[I[0]]
 
     def _batch_shape(self, shape=None, first_dim=False, B=None):
         if B is None:
@@ -300,25 +305,28 @@ def _init_dof_fields(self):
             ctrl_mode=gs.ti_int,
             hibernated=gs.ti_int,  # Flag for dofs that converge into a static state (hibernation)
         )
-        self.dofs_info = struct_dof_info.field(shape=self.n_dofs_, needs_grad=False, layout=ti.Layout.SOA)
+        dofs_info_shape = self._batch_shape(self.n_dofs_) if self._options.batch_dofs_info else self.n_dofs_
+        self.dofs_info = struct_dof_info.field(shape=dofs_info_shape, needs_grad=False, layout=ti.Layout.SOA)
         self.dofs_state = struct_dof_state.field(
             shape=self._batch_shape(self.n_dofs_), needs_grad=False, layout=ti.Layout.SOA
         )
 
         joints = self.joints
-        self._kernel_init_dof_fields(
-            dofs_motion_ang=np.concatenate([joint.dofs_motion_ang for joint in joints], dtype=gs.np_float),
-            dofs_motion_vel=np.concatenate([joint.dofs_motion_vel for joint in joints], dtype=gs.np_float),
-            dofs_limit=np.concatenate([joint.dofs_limit for joint in joints], dtype=gs.np_float),
-            dofs_invweight=np.concatenate([joint.dofs_invweight for joint in joints], dtype=gs.np_float),
-            dofs_stiffness=np.concatenate([joint.dofs_stiffness for joint in joints], dtype=gs.np_float),
-            dofs_sol_params=np.concatenate([joint.dofs_sol_params for joint in joints], dtype=gs.np_float),
-            dofs_damping=np.concatenate([joint.dofs_damping for joint in joints], dtype=gs.np_float),
-            dofs_armature=np.concatenate([joint.dofs_armature for joint in joints], dtype=gs.np_float),
-            dofs_kp=np.concatenate([joint.dofs_kp for joint in joints], dtype=gs.np_float),
-            dofs_kv=np.concatenate([joint.dofs_kv for joint in joints], dtype=gs.np_float),
-            dofs_force_range=np.concatenate([joint.dofs_force_range for joint in joints], dtype=gs.np_float),
-        )
+        is_nonempty = np.concatenate([joint.dofs_motion_ang for joint in joints], dtype=gs.np_float).shape[0] > 0
+        if is_nonempty:  # handle the case where there is a link with no dofs -- otherwise may cause invalid memory
+            self._kernel_init_dof_fields(
+                dofs_motion_ang=np.concatenate([joint.dofs_motion_ang for joint in joints], dtype=gs.np_float),
+                dofs_motion_vel=np.concatenate([joint.dofs_motion_vel for joint in joints], dtype=gs.np_float),
+                dofs_limit=np.concatenate([joint.dofs_limit for joint in joints], dtype=gs.np_float),
+                dofs_invweight=np.concatenate([joint.dofs_invweight for joint in joints], dtype=gs.np_float),
+                dofs_stiffness=np.concatenate([joint.dofs_stiffness for joint in joints], dtype=gs.np_float),
+                dofs_sol_params=np.concatenate([joint.dofs_sol_params for joint in joints], dtype=gs.np_float),
+                dofs_damping=np.concatenate([joint.dofs_damping for joint in joints], dtype=gs.np_float),
+                dofs_armature=np.concatenate([joint.dofs_armature for joint in joints], dtype=gs.np_float),
+                dofs_kp=np.concatenate([joint.dofs_kp for joint in joints], dtype=gs.np_float),
+                dofs_kv=np.concatenate([joint.dofs_kv for joint in joints], dtype=gs.np_float),
+                dofs_force_range=np.concatenate([joint.dofs_force_range for joint in joints], dtype=gs.np_float),
+            )
 
         # just in case
         self.dofs_state.force.fill(0)
@@ -339,26 +347,28 @@ def _kernel_init_dof_fields(
         dofs_force_range: ti.types.ndarray(),
     ):
         ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.PARTIAL)
-        for i, b in ti.ndrange(self.n_dofs, self._B):
+        for I in ti.grouped(self.dofs_info):
+            i = I[0]  # batching (if any) will be the second dim
+
             for j in ti.static(range(3)):
-                self.dofs_info[i].motion_ang[j] = dofs_motion_ang[i, j]
-                self.dofs_info[i].motion_vel[j] = dofs_motion_vel[i, j]
+                self.dofs_info[I].motion_ang[j] = dofs_motion_ang[i, j]
+                self.dofs_info[I].motion_vel[j] = dofs_motion_vel[i, j]
 
             for j in ti.static(range(2)):
-                self.dofs_info[i].limit[j] = dofs_limit[i, j]
-                self.dofs_info[i].force_range[j] = dofs_force_range[i, j]
+                self.dofs_info[I].limit[j] = dofs_limit[i, j]
+                self.dofs_info[I].force_range[j] = dofs_force_range[i, j]
 
             for j in ti.static(range(7)):
-                self.dofs_info[i].sol_params[j] = dofs_sol_params[i, j]
+                self.dofs_info[I].sol_params[j] = dofs_sol_params[i, j]
 
-            self.dofs_info[i].sol_params[0] = self._sol_contact_resolve_time
+            self.dofs_info[I].sol_params[0] = self._sol_contact_resolve_time
 
-            self.dofs_info[i].armature = dofs_armature[i]
-            self.dofs_info[i].invweight = dofs_invweight[i]
-            self.dofs_info[i].stiffness = dofs_stiffness[i]
-            self.dofs_info[i].damping = dofs_damping[i]
-            self.dofs_info[i].kp = dofs_kp[i]
-            self.dofs_info[i].kv = dofs_kv[i]
+            self.dofs_info[I].armature = dofs_armature[i]
+            self.dofs_info[I].invweight = dofs_invweight[i]
+            self.dofs_info[I].stiffness = dofs_stiffness[i]
+            self.dofs_info[I].damping = dofs_damping[i]
+            self.dofs_info[I].kp = dofs_kp[i]
+            self.dofs_info[I].kv = dofs_kv[i]
 
         ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.PARTIAL)
         for i, b in ti.ndrange(self.n_dofs, self._B):
@@ -442,7 +452,8 @@ def _init_link_fields(self):
             hibernated=gs.ti_int,
         )
 
-        self.links_info = struct_link_info.field(shape=self.n_links, needs_grad=False, layout=ti.Layout.SOA)
+        links_info_shape = self._batch_shape(self.n_links) if self._options.batch_links_info else self.n_links
+        self.links_info = struct_link_info.field(shape=links_info_shape, needs_grad=False, layout=ti.Layout.SOA)
         self.links_state = struct_link_state.field(
             shape=self._batch_shape(self.n_links), needs_grad=False, layout=ti.Layout.SOA
         )
@@ -501,38 +512,42 @@ def _kernel_init_link_fields(
         links_entity_idx: ti.types.ndarray(),
     ):
         ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.PARTIAL)
-        for i, b in ti.ndrange(self.n_links, self._B):
-            self.links_info[i].parent_idx = links_parent_idx[i]
-            self.links_info[i].root_idx = links_root_idx[i]
-            self.links_info[i].q_start = links_q_start[i]
-            self.links_info[i].dof_start = links_dof_start[i]
-            self.links_info[i].q_end = links_q_end[i]
-            self.links_info[i].dof_end = links_dof_end[i]
-            self.links_info[i].n_dofs = links_dof_end[i] - links_dof_start[i]
-            self.links_info[i].joint_type = links_joint_type[i]
-            self.links_info[i].invweight = links_invweight[i]
-            self.links_info[i].is_fixed = links_is_fixed[i]
-            self.links_info[i].entity_idx = links_entity_idx[i]
+        for I in ti.grouped(self.links_info):
+            i = I[0]
+
+            self.links_info[I].parent_idx = links_parent_idx[i]
+            self.links_info[I].root_idx = links_root_idx[i]
+            self.links_info[I].q_start = links_q_start[i]
+            self.links_info[I].dof_start = links_dof_start[i]
+            self.links_info[I].q_end = links_q_end[i]
+            self.links_info[I].dof_end = links_dof_end[i]
+            self.links_info[I].n_dofs = links_dof_end[i] - links_dof_start[i]
+            self.links_info[I].joint_type = links_joint_type[i]
+            self.links_info[I].invweight = links_invweight[i]
+            self.links_info[I].is_fixed = links_is_fixed[i]
+            self.links_info[I].entity_idx = links_entity_idx[i]
 
             for j in ti.static(range(4)):
-                self.links_info[i].quat[j] = links_quat[i, j]
-                self.links_info[i].joint_quat[j] = links_joint_quat[i, j]
-                self.links_info[i].inertial_quat[j] = links_inertial_quat[i, j]
+                self.links_info[I].quat[j] = links_quat[i, j]
+                self.links_info[I].joint_quat[j] = links_joint_quat[i, j]
+                self.links_info[I].inertial_quat[j] = links_inertial_quat[i, j]
 
             for j in ti.static(range(3)):
-                self.links_info[i].pos[j] = links_pos[i, j]
-                self.links_info[i].joint_pos[j] = links_joint_pos[i, j]
-                self.links_info[i].inertial_pos[j] = links_inertial_pos[i, j]
+                self.links_info[I].pos[j] = links_pos[i, j]
+                self.links_info[I].joint_pos[j] = links_joint_pos[i, j]
+                self.links_info[I].inertial_pos[j] = links_inertial_pos[i, j]
 
-            self.links_info[i].inertial_mass = links_inertial_mass[i]
+            self.links_info[I].inertial_mass = links_inertial_mass[i]
             for j1 in ti.static(range(3)):
                 for j2 in ti.static(range(3)):
-                    self.links_info[i].inertial_i[j1, j2] = links_inertial_i[i, j1, j2]
+                    self.links_info[I].inertial_i[j1, j2] = links_inertial_i[i, j1, j2]
 
         ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.PARTIAL)
         for i, b in ti.ndrange(self.n_links, self._B):
+            I = [i, b] if ti.static(self._options.batch_links_info) else i
+
             # Update state for root fixed link. Their state will not be updated in forward kinematics later but can be manually changed by user.
-            if self.links_info[i].parent_idx == -1 and self.links_info[i].is_fixed:
+            if self.links_info[I].parent_idx == -1 and self.links_info[I].is_fixed:
                 for j in ti.static(range(4)):
                     self.links_state[i, b].quat[j] = links_quat[i, j]
 
@@ -836,10 +851,11 @@ def _kernel_adjust_link_inertia(
         link_idx: ti.i32,
         ratio: ti.f32,
     ):
-        self.links_info[link_idx].invweight /= ratio
-        self.links_info[link_idx].inertial_mass *= ratio
-        for j1, j2 in ti.ndrange(3, 3):
-            self.links_info[link_idx].inertial_i[j1, j2] *= ratio
+        for I_l in ti.grouped(self.links_info):
+            self.links_info[I_l].invweight /= ratio
+            self.links_info[I_l].inertial_mass *= ratio
+            for j1, j2 in ti.ndrange(3, 3):
+                self.links_info[I_l].inertial_i[j1, j2] *= ratio
 
     def _init_vgeom_fields(self):
         struct_vgeom_info = ti.types.struct(
@@ -971,8 +987,13 @@ def _kernel_init_entity_fields(
 
             self.entities_info[i].gravity_compensation = entities_gravity_compensation[i]
 
-            for i_d in range(entities_dof_start[i], entities_dof_end[i]):
-                self.dofs_info[i_d].dof_start = entities_dof_start[i]
+            if ti.static(self._options.batch_dofs_info):
+                for i_b in range(self._B):
+                    for i_d in range(entities_dof_start[i], entities_dof_end[i]):
+                        self.dofs_info[i_d, i_b].dof_start = entities_dof_start[i]
+            else:
+                for i_d in range(entities_dof_start[i], entities_dof_end[i]):
+                    self.dofs_info[i_d].dof_start = entities_dof_start[i]
 
         if ti.static(self._use_hibernation):
             ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.PARTIAL)
@@ -1053,7 +1074,8 @@ def _func_compute_mass_matrix(self):
                     i_e = self.awake_entities[i_e_, i_b]
                     for i in range(self.entities_info[i_e].n_links):
                         i_l = self.entities_info[i_e].link_end - 1 - i
-                        i_p = self.links_info[i_l].parent_idx
+                        I_l = [i_l, i_b] if ti.static(self._options.batch_links_info) else i_l
+                        i_p = self.links_info[I_l].parent_idx
 
                         if i_p != -1:
                             self.links_state[i_p, i_b].crb_inertial = (
@@ -1075,7 +1097,8 @@ def _func_compute_mass_matrix(self):
             for i_b in range(self._B):
                 for i_l_ in range(self.n_awake_links[i_b]):
                     i_l = self.awake_links[i_l_, i_b]
-                    l_info = self.links_info[i_l]
+                    I_l = [i_l, i_b] if ti.static(self._options.batch_links_info) else i_l
+                    l_info = self.links_info[I_l]
                     for i_d in range(l_info.dof_start, l_info.dof_end):
                         self.dofs_state[i_d, i_b].f_ang, self.dofs_state[i_d, i_b].f_vel = gu.inertial_mul(
                             self.links_state[i_l, i_b].crb_pos,
@@ -1098,10 +1121,11 @@ def _func_compute_mass_matrix(self):
                             ) * self.mass_parent_mask[i_d, j_d]
 
                     for i_d in range(e_info.dof_start, e_info.dof_end):
+                        I_d = [i_d, i_b] if ti.static(self._options.batch_dofs_info) else i_d
                         self.mass_mat[i_d, i_d, i_b] = (
                             self.mass_mat[i_d, i_d, i_b]
-                            + self.dofs_info[i_d].armature
-                            + self.dofs_info[i_d].damping * self._substep_dt
+                            + self.dofs_info[I_d].armature
+                            + self.dofs_info[I_d].damping * self._substep_dt
                         )
                         for j_d in range(i_d + 1, e_info.dof_end):
                             self.mass_mat[i_d, j_d, i_b] = self.mass_mat[j_d, i_d, i_b]
@@ -1111,18 +1135,19 @@ def _func_compute_mass_matrix(self):
                 # qDeriv += d qfrc_actuator / d qvel
                 ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.ALL)
                 for i_d, i_b in ti.ndrange(self.n_dofs, self._B):
+                    I_d = [i_d, i_b] if ti.static(self._options.batch_dofs_info) else i_d
 
                     if self.dofs_state[i_d, i_b].ctrl_mode == gs.CTRL_MODE.FORCE:
                         pass
 
                     elif self.dofs_state[i_d, i_b].ctrl_mode == gs.CTRL_MODE.VELOCITY:
                         self.mass_mat[i_d, i_d, i_b] = (
-                            self.mass_mat[i_d, i_d, i_b] + self.dofs_info[i_d].kv * self._substep_dt
+                            self.mass_mat[i_d, i_d, i_b] + self.dofs_info[I_d].kv * self._substep_dt
                         )
 
                     elif self.dofs_state[i_d, i_b].ctrl_mode == gs.CTRL_MODE.POSITION:
                         self.mass_mat[i_d, i_d, i_b] = (
-                            self.mass_mat[i_d, i_d, i_b] + self.dofs_info[i_d].kv * self._substep_dt
+                            self.mass_mat[i_d, i_d, i_b] + self.dofs_info[I_d].kv * self._substep_dt
                         )
 
         else:
@@ -1139,7 +1164,8 @@ def _func_compute_mass_matrix(self):
             for i_e, i_b in ti.ndrange(self.n_entities, self._B):
                 for i in range(self.entities_info[i_e].n_links):
                     i_l = self.entities_info[i_e].link_end - 1 - i
-                    i_p = self.links_info[i_l].parent_idx
+                    I_l = [i_l, i_b] if ti.static(self._options.batch_links_info) else i_l
+                    i_p = self.links_info[I_l].parent_idx
 
                     if i_p != -1:
                         self.links_state[i_p, i_b].crb_inertial = (
@@ -1160,7 +1186,8 @@ def _func_compute_mass_matrix(self):
             ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.ALL)
             for i_b in range(self._B):
                 for i_l in range(self.n_links):
-                    l_info = self.links_info[i_l]
+                    I_l = [i_l, i_b] if ti.static(self._options.batch_links_info) else i_l
+                    l_info = self.links_info[I_l]
                     for i_d in range(l_info.dof_start, l_info.dof_end):
                         self.dofs_state[i_d, i_b].f_ang, self.dofs_state[i_d, i_b].f_vel = gu.inertial_mul(
                             self.links_state[i_l, i_b].crb_pos,
@@ -1181,10 +1208,11 @@ def _func_compute_mass_matrix(self):
                         ) * self.mass_parent_mask[i_d, j_d]
 
                 for i_d in range(e_info.dof_start, e_info.dof_end):
+                    I_d = [i_d, i_b] if ti.static(self._options.batch_dofs_info) else i_d
                     self.mass_mat[i_d, i_d, i_b] = (
                         self.mass_mat[i_d, i_d, i_b]
-                        + self.dofs_info[i_d].armature
-                        + self.dofs_info[i_d].damping * self._substep_dt
+                        + self.dofs_info[I_d].armature
+                        + self.dofs_info[I_d].damping * self._substep_dt
                     )
                     for j_d in range(i_d + 1, e_info.dof_end):
                         self.mass_mat[i_d, j_d, i_b] = self.mass_mat[j_d, i_d, i_b]
@@ -1194,18 +1222,19 @@ def _func_compute_mass_matrix(self):
                 # qDeriv += d qfrc_actuator / d qvel
                 ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.ALL)
                 for i_d, i_b in ti.ndrange(self.n_dofs, self._B):
+                    I_d = [i_d, i_b] if ti.static(self._options.batch_dofs_info) else i_d
 
                     if self.dofs_state[i_d, i_b].ctrl_mode == gs.CTRL_MODE.FORCE:
                         pass
 
                     elif self.dofs_state[i_d, i_b].ctrl_mode == gs.CTRL_MODE.VELOCITY:
                         self.mass_mat[i_d, i_d, i_b] = (
-                            self.mass_mat[i_d, i_d, i_b] + self.dofs_info[i_d].kv * self._substep_dt
+                            self.mass_mat[i_d, i_d, i_b] + self.dofs_info[I_d].kv * self._substep_dt
                         )
 
                     elif self.dofs_state[i_d, i_b].ctrl_mode == gs.CTRL_MODE.POSITION:
                         self.mass_mat[i_d, i_d, i_b] = (
-                            self.mass_mat[i_d, i_d, i_b] + self.dofs_info[i_d].kv * self._substep_dt
+                            self.mass_mat[i_d, i_d, i_b] + self.dofs_info[I_d].kv * self._substep_dt
                         )
 
     @ti.func
@@ -1353,15 +1382,16 @@ def _func_implicit_damping(self):
         # TODO: hibernate
         ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.ALL)
         for i_d, i_b in ti.ndrange(self.n_dofs, self._B):
+            I_d = [i_d, i_b] if ti.static(self._options.batch_dofs_info) else i_d
 
             if self.dofs_state[i_d, i_b].ctrl_mode == gs.CTRL_MODE.FORCE:
                 pass
 
             elif self.dofs_state[i_d, i_b].ctrl_mode == gs.CTRL_MODE.VELOCITY:
-                self.mass_mat[i_d, i_d, i_b] = self.mass_mat[i_d, i_d, i_b] + self.dofs_info[i_d].kv * self._substep_dt
+                self.mass_mat[i_d, i_d, i_b] = self.mass_mat[i_d, i_d, i_b] + self.dofs_info[I_d].kv * self._substep_dt
 
             elif self.dofs_state[i_d, i_b].ctrl_mode == gs.CTRL_MODE.POSITION:
-                self.mass_mat[i_d, i_d, i_b] = self.mass_mat[i_d, i_d, i_b] + self.dofs_info[i_d].kv * self._substep_dt
+                self.mass_mat[i_d, i_d, i_b] = self.mass_mat[i_d, i_d, i_b] + self.dofs_info[I_d].kv * self._substep_dt
 
             self.dofs_state[i_d, i_b].force += self.dofs_state[i_d, i_b].qf_constraint
 
@@ -1491,9 +1521,10 @@ def _func_COM_links(self):
             for i_b in range(self._B):
                 for i_l_ in range(self.n_awake_links[i_b]):
                     i_l = self.awake_links[i_l_, i_b]
+                    I_l = [i_l, i_b] if ti.static(self._options.batch_links_info) else i_l
 
                     l = self.links_state[i_l, i_b]
-                    l_info = self.links_info[i_l]
+                    l_info = self.links_info[I_l]
                     mass = l_info.inertial_mass + l.mass_shift
                     self.links_state[i_l, i_b].i_pos, self.links_state[i_l, i_b].i_quat = (
                         gu.ti_transform_pos_quat_by_trans_quat(
@@ -1501,7 +1532,7 @@ def _func_COM_links(self):
                         )
                     )
 
-                    i_r = self.links_info[i_l].root_idx
+                    i_r = self.links_info[I_l].root_idx
                     ti.atomic_add(self.links_state[i_r, i_b].mass_sum, mass)
 
                     COM = mass * self.links_state[i_l, i_b].i_pos
@@ -1511,8 +1542,9 @@ def _func_COM_links(self):
             for i_b in range(self._B):
                 for i_l_ in range(self.n_awake_links[i_b]):
                     i_l = self.awake_links[i_l_, i_b]
+                    I_l = [i_l, i_b] if ti.static(self._options.batch_links_info) else i_l
 
-                    i_r = self.links_info[i_l].root_idx
+                    i_r = self.links_info[I_l].root_idx
                     if i_l == i_r:
                         self.links_state[i_l, i_b].root_COM = (
                             self.links_state[i_l, i_b].root_COM / self.links_state[i_l, i_b].mass_sum
@@ -1522,19 +1554,21 @@ def _func_COM_links(self):
             for i_b in range(self._B):
                 for i_l_ in range(self.n_awake_links[i_b]):
                     i_l = self.awake_links[i_l_, i_b]
+                    I_l = [i_l, i_b] if ti.static(self._options.batch_links_info) else i_l
 
-                    i_r = self.links_info[i_l].root_idx
+                    i_r = self.links_info[I_l].root_idx
                     self.links_state[i_l, i_b].root_COM = self.links_state[i_r, i_b].root_COM
 
             ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.ALL)
             for i_b in range(self._B):
                 for i_l_ in range(self.n_awake_links[i_b]):
                     i_l = self.awake_links[i_l_, i_b]
+                    I_l = [i_l, i_b] if ti.static(self._options.batch_links_info) else i_l
 
                     l = self.links_state[i_l, i_b]
-                    l_info = self.links_info[i_l]
+                    l_info = self.links_info[I_l]
 
-                    i_r = self.links_info[i_l].root_idx
+                    i_r = self.links_info[I_l].root_idx
                     self.links_state[i_l, i_b].COM = self.links_state[i_r, i_b].root_COM
                     self.links_state[i_l, i_b].i_pos = self.links_state[i_l, i_b].i_pos - self.links_state[i_l, i_b].COM
 
@@ -1553,8 +1587,9 @@ def _func_COM_links(self):
             for i_b in range(self._B):
                 for i_l_ in range(self.n_awake_links[i_b]):
                     i_l = self.awake_links[i_l_, i_b]
+                    I_l = [i_l, i_b] if ti.static(self._options.batch_links_info) else i_l
 
-                    l_info = self.links_info[i_l]
+                    l_info = self.links_info[I_l]
                     i_p = l_info.parent_idx
 
                     p_pos = ti.Vector.zero(gs.ti_float, 3)
@@ -1586,14 +1621,16 @@ def _func_COM_links(self):
             for i_b in range(self._B):
                 for i_l_ in range(self.n_awake_links[i_b]):
                     i_l = self.awake_links[i_l_, i_b]
+                    I_l = [i_l, i_b] if ti.static(self._options.batch_links_info) else i_l
 
-                    l_info = self.links_info[i_l]
+                    l_info = self.links_info[I_l]
 
                     if l_info.joint_type == gs.JOINT_TYPE.FREE:
                         for i_d in range(l_info.dof_start, l_info.dof_end):
-                            self.dofs_state[i_d, i_b].cdof_vel = self.dofs_info[i_d].motion_vel
+                            I_d = [i_d, i_b] if ti.static(self._options.batch_dofs_info) else i_d
+                            self.dofs_state[i_d, i_b].cdof_vel = self.dofs_info[I_d].motion_vel
                             self.dofs_state[i_d, i_b].cdof_ang = gu.ti_transform_by_quat(
-                                self.dofs_info[i_d].motion_ang, self.links_state[i_l, i_b].j_quat
+                                self.dofs_info[I_d].motion_ang, self.links_state[i_l, i_b].j_quat
                             )
 
                             offset_pos = self.links_state[i_l, i_b].COM - self.links_state[i_l, i_b].j_pos
@@ -1617,8 +1654,9 @@ def _func_COM_links(self):
                         pass
                     else:
                         for i_d in range(l_info.dof_start, l_info.dof_end):
-                            motion_vel = self.dofs_info[i_d].motion_vel
-                            motion_ang = self.dofs_info[i_d].motion_ang
+                            I_d = [i_d, i_b] if ti.static(self._options.batch_dofs_info) else i_d
+                            motion_vel = self.dofs_info[I_d].motion_vel
+                            motion_ang = self.dofs_info[I_d].motion_ang
 
                             self.dofs_state[i_d, i_b].cdof_ang = gu.ti_transform_by_quat(
                                 motion_ang, self.links_state[i_l, i_b].j_quat
@@ -1654,9 +1692,10 @@ def _func_COM_links(self):
             ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.ALL)
             for i_b in range(self._B):
                 for i_l in range(self.n_links):
+                    I_l = [i_l, i_b] if ti.static(self._options.batch_links_info) else i_l
 
                     l = self.links_state[i_l, i_b]
-                    l_info = self.links_info[i_l]
+                    l_info = self.links_info[I_l]
                     mass = l_info.inertial_mass + l.mass_shift
                     self.links_state[i_l, i_b].i_pos, self.links_state[i_l, i_b].i_quat = (
                         gu.ti_transform_pos_quat_by_trans_quat(
@@ -1664,7 +1703,7 @@ def _func_COM_links(self):
                         )
                     )
 
-                    i_r = self.links_info[i_l].root_idx
+                    i_r = self.links_info[I_l].root_idx
                     ti.atomic_add(self.links_state[i_r, i_b].mass_sum, mass)
 
                     COM = mass * self.links_state[i_l, i_b].i_pos
@@ -1673,8 +1712,9 @@ def _func_COM_links(self):
             ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.ALL)
             for i_b in range(self._B):
                 for i_l in range(self.n_links):
+                    I_l = [i_l, i_b] if ti.static(self._options.batch_links_info) else i_l
 
-                    i_r = self.links_info[i_l].root_idx
+                    i_r = self.links_info[I_l].root_idx
                     if i_l == i_r:
                         self.links_state[i_l, i_b].root_COM = (
                             self.links_state[i_l, i_b].root_COM / self.links_state[i_l, i_b].mass_sum
@@ -1683,18 +1723,20 @@ def _func_COM_links(self):
             ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.ALL)
             for i_b in range(self._B):
                 for i_l in range(self.n_links):
+                    I_l = [i_l, i_b] if ti.static(self._options.batch_links_info) else i_l
 
-                    i_r = self.links_info[i_l].root_idx
+                    i_r = self.links_info[I_l].root_idx
                     self.links_state[i_l, i_b].root_COM = self.links_state[i_r, i_b].root_COM
 
             ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.ALL)
             for i_b in range(self._B):
                 for i_l in range(self.n_links):
+                    I_l = [i_l, i_b] if ti.static(self._options.batch_links_info) else i_l
 
                     l = self.links_state[i_l, i_b]
-                    l_info = self.links_info[i_l]
+                    l_info = self.links_info[I_l]
 
-                    i_r = self.links_info[i_l].root_idx
+                    i_r = self.links_info[I_l].root_idx
                     self.links_state[i_l, i_b].COM = self.links_state[i_r, i_b].root_COM
                     self.links_state[i_l, i_b].i_pos = self.links_state[i_l, i_b].i_pos - self.links_state[i_l, i_b].COM
 
@@ -1712,8 +1754,9 @@ def _func_COM_links(self):
             ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.ALL)
             for i_b in range(self._B):
                 for i_l in range(self.n_links):
+                    I_l = [i_l, i_b] if ti.static(self._options.batch_links_info) else i_l
 
-                    l_info = self.links_info[i_l]
+                    l_info = self.links_info[I_l]
                     i_p = l_info.parent_idx
 
                     p_pos = ti.Vector.zero(gs.ti_float, 3)
@@ -1744,14 +1787,16 @@ def _func_COM_links(self):
             ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.ALL)
             for i_b in range(self._B):
                 for i_l in range(self.n_links):
+                    I_l = [i_l, i_b] if ti.static(self._options.batch_links_info) else i_l
 
-                    l_info = self.links_info[i_l]
+                    l_info = self.links_info[I_l]
 
                     if l_info.joint_type == gs.JOINT_TYPE.FREE:
                         for i_d in range(l_info.dof_start, l_info.dof_end):
-                            self.dofs_state[i_d, i_b].cdof_vel = self.dofs_info[i_d].motion_vel
+                            I_d = [i_d, i_b] if ti.static(self._options.batch_dofs_info) else i_d
+                            self.dofs_state[i_d, i_b].cdof_vel = self.dofs_info[I_d].motion_vel
                             self.dofs_state[i_d, i_b].cdof_ang = gu.ti_transform_by_quat(
-                                self.dofs_info[i_d].motion_ang, self.links_state[i_l, i_b].j_quat
+                                self.dofs_info[I_d].motion_ang, self.links_state[i_l, i_b].j_quat
                             )
 
                             offset_pos = self.links_state[i_l, i_b].COM - self.links_state[i_l, i_b].j_pos
@@ -1776,8 +1821,9 @@ def _func_COM_links(self):
 
                     else:
                         for i_d in range(l_info.dof_start, l_info.dof_end):
-                            motion_vel = self.dofs_info[i_d].motion_vel
-                            motion_ang = self.dofs_info[i_d].motion_ang
+                            I_d = [i_d, i_b] if ti.static(self._options.batch_dofs_info) else i_d
+                            motion_vel = self.dofs_info[I_d].motion_vel
+                            motion_ang = self.dofs_info[I_d].motion_ang
 
                             self.dofs_state[i_d, i_b].cdof_ang = gu.ti_transform_by_quat(
                                 motion_ang, self.links_state[i_l, i_b].j_quat
@@ -1836,7 +1882,8 @@ def _func_COM_cd(self):
 
                 e_info = self.entities_info[i_e]
                 for i_l in range(e_info.link_start, e_info.link_end):
-                    l_info = self.links_info[i_l]
+                    I_l = [i_l, i_b] if ti.static(self._options.batch_links_info) else i_l
+                    l_info = self.links_info[I_l]
                     i_p = l_info.parent_idx
 
                     cd_vel = ti.Vector.zero(gs.ti_float, 3)
@@ -1861,8 +1908,9 @@ def _func_COM_cdofd(self):
             for i_b in range(self._B):
                 for i_l_ in range(self.n_awake_links[i_b]):
                     i_l = self.awake_links[i_l_, i_b]
+                    I_l = [i_l, i_b] if ti.static(self._options.batch_links_info) else i_l
 
-                    l_info = self.links_info[i_l]
+                    l_info = self.links_info[I_l]
 
                     if l_info.joint_type == gs.JOINT_TYPE.FREE:
                         cd_ang = ti.Vector.zero(gs.ti_float, 3)
@@ -1904,8 +1952,9 @@ def _func_COM_cdofd(self):
             ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.ALL)
             for i_b in range(self._B):
                 for i_l in range(self.n_links):
+                    I_l = [i_l, i_b] if ti.static(self._options.batch_links_info) else i_l
 
-                    l_info = self.links_info[i_l]
+                    l_info = self.links_info[I_l]
 
                     if l_info.joint_type == gs.JOINT_TYPE.FREE:
                         cd_ang = ti.Vector.zero(gs.ti_float, 3)
@@ -1967,7 +2016,8 @@ def _func_forward_kinematics(self):
     def _func_forward_kinematics_entity(self, i_e, i_b):
         # calculate_j
         for i_l in range(self.entities_info[i_e].link_start, self.entities_info[i_e].link_end):
-            l_info = self.links_info[i_l]
+            I_l = [i_l, i_b] if ti.static(self._options.batch_links_info) else i_l
+            l_info = self.links_info[I_l]
 
             if l_info.joint_type == gs.JOINT_TYPE.FREE:
 
@@ -1991,7 +2041,8 @@ def _func_forward_kinematics_entity(self, i_e, i_b):
 
             else:
                 self.dofs_state[l_info.dof_start, i_b].pos = self.qpos[l_info.q_start, i_b]
-                dof_info = self.dofs_info[l_info.dof_start]
+                I_dof_start = [l_info.dof_start, i_b] if ti.static(self._options.batch_dofs_info) else l_info.dof_start
+                dof_info = self.dofs_info[I_dof_start]
                 self.links_state[i_l, i_b].j_pos = dof_info.motion_vel * self.qpos[l_info.q_start, i_b]
                 self.links_state[i_l, i_b].j_quat = gu.ti_rotvec_to_quat(
                     dof_info.motion_ang * self.qpos[l_info.q_start, i_b]
@@ -2000,7 +2051,8 @@ def _func_forward_kinematics_entity(self, i_e, i_b):
                 self.links_state[i_l, i_b].j_vel = dof_info.motion_vel * self.dofs_state[l_info.dof_start, i_b].vel
 
                 for i_d in range(l_info.dof_start + 1, l_info.dof_end):
-                    dof_info = self.dofs_info[i_d]
+                    I_d = [i_d, i_b] if ti.static(self._options.batch_dofs_info) else i_d
+                    dof_info = self.dofs_info[I_d]
                     qi = l_info.q_start + i_d - l_info.dof_start
                     self.dofs_state[i_d, i_b].pos = self.qpos[qi, i_b]
                     ji_pos = dof_info.motion_vel * self.qpos[qi, i_b]
@@ -2041,7 +2093,8 @@ def _func_forward_kinematics_entity(self, i_e, i_b):
         # joint_to_world
         for i_l in range(self.entities_info[i_e].link_start, self.entities_info[i_e].link_end):
             l_state = self.links_state[i_l, i_b]
-            l_info = self.links_info[i_l]
+            I_l = [i_l, i_b] if ti.static(self._options.batch_links_info) else i_l
+            l_info = self.links_info[I_l]
             i_p = l_info.parent_idx
 
             if i_p == -1:  # root link
@@ -2318,15 +2371,17 @@ def _func_apply_external_torque_link_frame(self, torque, link_idx, batch_idx):
 
     @ti.func
     def _func_apply_external_force_link_inertial_frame(self, pos, force, link_idx, batch_idx):
+        link_I = [link_idx, batch_idx] if ti.static(self._options.batch_links_info) else link_idx
         pos = gu.ti_transform_by_trans_quat(
-            pos, self.links_info[link_idx].inertial_pos, self.links_info[link_idx].inertial_quat
+            pos, self.links_info[link_I].inertial_pos, self.links_info[link_I].inertial_quat
         )
-        force = gu.ti_transform_by_quat(force, self.links_info[link_idx].inertial_quat)
+        force = gu.ti_transform_by_quat(force, self.links_info[link_I].inertial_quat)
         self._func_apply_external_force_link_frame(pos, force, link_idx, batch_idx)
 
     @ti.func
     def _func_apply_external_torque_link_inertial_frame(self, torque, link_idx, batch_idx):
-        torque = gu.ti_transform_by_quat(torque, self.links_info[link_idx].inertial_quat)
+        link_I = [link_idx, batch_idx] if ti.static(self._options.batch_links_info) else link_idx
+        torque = gu.ti_transform_by_quat(torque, self.links_info[link_I].inertial_quat)
         self._func_apply_external_torque_link_frame(torque, link_idx, batch_idx)
 
     @ti.func
@@ -2352,27 +2407,29 @@ def _func_torque_and_passive_force(self):
             wakeup = False
             for i_l in range(self.entities_info[i_e].link_start, self.entities_info[i_e].link_end):
                 force = gs.ti_float(0.0)
-                l_info = self.links_info[i_l]
+                I_l = [i_l, i_b] if ti.static(self._options.batch_links_info) else i_l
+                l_info = self.links_info[I_l]
                 for i_d in range(l_info.dof_start, l_info.dof_end):
+                    I_d = [i_d, i_b] if ti.static(self._options.batch_dofs_info) else i_d
                     if self.dofs_state[i_d, i_b].ctrl_mode == gs.CTRL_MODE.FORCE:
                         force = self.dofs_state[i_d, i_b].ctrl_force
                     elif self.dofs_state[i_d, i_b].ctrl_mode == gs.CTRL_MODE.VELOCITY:
-                        force = self.dofs_info[i_d].kv * (
+                        force = self.dofs_info[I_d].kv * (
                             self.dofs_state[i_d, i_b].ctrl_vel - self.dofs_state[i_d, i_b].vel
                         )
                     elif self.dofs_state[i_d, i_b].ctrl_mode == gs.CTRL_MODE.POSITION and not (
                         l_info.joint_type == gs.JOINT_TYPE.FREE and i_d >= l_info.dof_start + 3
                     ):
                         force = (
-                            self.dofs_info[i_d].kp
+                            self.dofs_info[I_d].kp
                             * (self.dofs_state[i_d, i_b].ctrl_pos - self.dofs_state[i_d, i_b].pos)
-                            - self.dofs_info[i_d].kv * self.dofs_state[i_d, i_b].vel
+                            - self.dofs_info[I_d].kv * self.dofs_state[i_d, i_b].vel
                         )
 
                     self.dofs_state[i_d, i_b].qf_applied = ti.math.clamp(
                         force,
-                        self.dofs_info[i_d].force_range[0],
-                        self.dofs_info[i_d].force_range[1],
+                        self.dofs_info[I_d].force_range[0],
+                        self.dofs_info[I_d].force_range[1],
                     )
 
                     if ti.abs(force) > gs.EPS:
@@ -2410,13 +2467,14 @@ def _func_torque_and_passive_force(self):
                     rotvec = gu.ti_quat_to_rotvec(q_diff)
 
                     for i_d in range(l_info.dof_start + 3, l_info.dof_end):
+                        I_d = [i_d, i_b] if ti.static(self._options.batch_dofs_info) else i_d
                         force = (
-                            self.dofs_info[i_d].kp * rotvec[i_d - l_info.dof_start - 3]
-                            - self.dofs_info[i_d].kv * self.dofs_state[i_d, i_b].vel
+                            self.dofs_info[I_d].kp * rotvec[i_d - l_info.dof_start - 3]
+                            - self.dofs_info[I_d].kv * self.dofs_state[i_d, i_b].vel
                         )
 
                         self.dofs_state[i_d, i_b].qf_applied = ti.math.clamp(
-                            force, self.dofs_info[i_d].force_range[0], self.dofs_info[i_d].force_range[1]
+                            force, self.dofs_info[I_d].force_range[0], self.dofs_info[I_d].force_range[1]
                         )
                         if ti.abs(force) > gs.EPS:
                             wakeup = True
@@ -2430,7 +2488,8 @@ def _func_torque_and_passive_force(self):
             for i_b in range(self._B):
                 for i_l_ in range(self.n_awake_links[i_b]):
                     i_l = self.awake_links[i_l_, i_b]
-                    l_info = self.links_info[i_l]
+                    I_l = [i_l, i_b] if ti.static(self._options.batch_links_info) else i_l
+                    l_info = self.links_info[I_l]
                     if (
                         l_info.joint_type == gs.JOINT_TYPE.REVOLUTE
                         or l_info.joint_type == gs.JOINT_TYPE.PRISMATIC
@@ -2442,20 +2501,25 @@ def _func_torque_and_passive_force(self):
                         q_end = l_info.q_end
 
                         for j_d in range(q_end - q_start):
+                            I_d = (
+                                [dof_start + j_d, i_b] if ti.static(self._options.batch_dofs_info) else dof_start + j_d
+                            )
                             self.dofs_state[dof_start + j_d, i_b].qf_passive = (
-                                -self.qpos[q_start + j_d, i_b] * self.dofs_info[dof_start + j_d].stiffness
+                                -self.qpos[q_start + j_d, i_b] * self.dofs_info[I_d].stiffness
                             )
 
             ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.ALL)
             for i_b in range(self._B):
                 for i_d_ in range(self.n_awake_dofs[i_b]):
                     i_d = self.awake_dofs[i_d_, i_b]
+                    I_d = [i_d, i_b] if ti.static(self._options.batch_dofs_info) else i_d
 
-                    self.dofs_state[i_d, i_b].qf_passive += -self.dofs_info[i_d].damping * self.dofs_state[i_d, i_b].vel
+                    self.dofs_state[i_d, i_b].qf_passive += -self.dofs_info[I_d].damping * self.dofs_state[i_d, i_b].vel
         else:
             ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.ALL)
             for i_l, i_b in ti.ndrange(self.n_links, self._B):
-                l_info = self.links_info[i_l]
+                I_l = [i_l, i_b] if ti.static(self._options.batch_links_info) else i_l
+                l_info = self.links_info[I_l]
                 if (
                     l_info.joint_type == gs.JOINT_TYPE.REVOLUTE
                     or l_info.joint_type == gs.JOINT_TYPE.PRISMATIC
@@ -2467,13 +2531,15 @@ def _func_torque_and_passive_force(self):
                     q_end = l_info.q_end
 
                     for j_d in range(q_end - q_start):
+                        I_d = [dof_start + j_d, i_b] if ti.static(self._options.batch_dofs_info) else dof_start + j_d
                         self.dofs_state[dof_start + j_d, i_b].qf_passive = (
-                            -self.qpos[q_start + j_d, i_b] * self.dofs_info[dof_start + j_d].stiffness
+                            -self.qpos[q_start + j_d, i_b] * self.dofs_info[I_d].stiffness
                         )
 
             ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.ALL)
             for i_d, i_b in ti.ndrange(self.n_dofs, self._B):
-                self.dofs_state[i_d, i_b].qf_passive += -self.dofs_info[i_d].damping * self.dofs_state[i_d, i_b].vel
+                I_d = [i_d, i_b] if ti.static(self._options.batch_dofs_info) else i_d
+                self.dofs_state[i_d, i_b].qf_passive += -self.dofs_info[I_d].damping * self.dofs_state[i_d, i_b].vel
 
     @ti.func
     def _func_system_update_acc(self):
@@ -2484,8 +2550,9 @@ def _func_system_update_acc(self):
                     i_e = self.awake_entities[i_e_, i_b]
                     e_info = self.entities_info[i_e]
                     for i_l in range(e_info.link_start, e_info.link_end):
+                        I_l = [i_l, i_b] if ti.static(self._options.batch_links_info) else i_l
 
-                        i_p = self.links_info[i_l].parent_idx
+                        i_p = self.links_info[I_l].parent_idx
                         if i_p == -1:
                             self.links_state[i_l, i_b].cdd_vel = -self._gravity[None] * (
                                 1 - e_info.gravity_compensation
@@ -2499,7 +2566,7 @@ def _func_system_update_acc(self):
                         map_sum_vel = ti.Vector.zero(gs.ti_float, 3)
                         map_sum_ang = ti.Vector.zero(gs.ti_float, 3)
 
-                        for i_d in range(self.links_info[i_l].dof_start, self.links_info[i_l].dof_end):
+                        for i_d in range(self.links_info[I_l].dof_start, self.links_info[I_l].dof_end):
                             map_sum_vel = (
                                 map_sum_vel + self.dofs_state[i_d, i_b].cdofd_vel * self.dofs_state[i_d, i_b].vel
                             )
@@ -2514,8 +2581,9 @@ def _func_system_update_acc(self):
             for i_e, i_b in ti.ndrange(self.n_entities, self._B):
                 e_info = self.entities_info[i_e]
                 for i_l in range(e_info.link_start, e_info.link_end):
+                    I_l = [i_l, i_b] if ti.static(self._options.batch_links_info) else i_l
 
-                    i_p = self.links_info[i_l].parent_idx
+                    i_p = self.links_info[I_l].parent_idx
                     if i_p == -1:
                         self.links_state[i_l, i_b].cdd_vel = -self._gravity[None] * (1 - e_info.gravity_compensation)
                         self.links_state[i_l, i_b].cdd_ang = ti.Vector.zero(gs.ti_float, 3)
@@ -2527,7 +2595,7 @@ def _func_system_update_acc(self):
                     map_sum_vel = ti.Vector.zero(gs.ti_float, 3)
                     map_sum_ang = ti.Vector.zero(gs.ti_float, 3)
 
-                    for i_d in range(self.links_info[i_l].dof_start, self.links_info[i_l].dof_end):
+                    for i_d in range(self.links_info[I_l].dof_start, self.links_info[I_l].dof_end):
                         map_sum_vel = map_sum_vel + self.dofs_state[i_d, i_b].cdofd_vel * self.dofs_state[i_d, i_b].vel
                         map_sum_ang = map_sum_ang + self.dofs_state[i_d, i_b].cdofd_ang * self.dofs_state[i_d, i_b].vel
 
@@ -2596,7 +2664,8 @@ def _func_inverse_link_force(self):
                     e_info = self.entities_info[i_e]
                     for i in range(e_info.n_links):
                         i_l = e_info.link_end - 1 - i
-                        i_p = self.links_info[i_l].parent_idx
+                        I_l = [i_l, i_b] if ti.static(self._options.batch_links_info) else i_l
+                        i_p = self.links_info[I_l].parent_idx
                         if i_p != -1:
                             self.links_state[i_p, i_b].cfrc_flat_vel = (
                                 self.links_state[i_p, i_b].cfrc_flat_vel + self.links_state[i_l, i_b].cfrc_flat_vel
@@ -2611,7 +2680,8 @@ def _func_inverse_link_force(self):
                 e_info = self.entities_info[i_e]
                 for i in range(e_info.n_links):
                     i_l = e_info.link_end - 1 - i
-                    i_p = self.links_info[i_l].parent_idx
+                    I_l = [i_l, i_b] if ti.static(self._options.batch_links_info) else i_l
+                    i_p = self.links_info[I_l].parent_idx
                     if i_p != -1:
                         self.links_state[i_p, i_b].cfrc_flat_vel = (
                             self.links_state[i_p, i_b].cfrc_flat_vel + self.links_state[i_l, i_b].cfrc_flat_vel
@@ -2628,15 +2698,16 @@ def _func_actuation(self):
             ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.ALL)
             for i_l, i_b in ti.ndrange(self.n_links, self._B):
                 joint_type = self.links_info[i_l].joint_type
-                q_start = self.links_info[i_l].q_start
+                I_l = [i_l, i_b] if ti.static(self._options.batch_links_info) else i_l
+                q_start = self.links_info[I_l].q_start
 
                 if joint_type == gs.JOINT_TYPE.REVOLUTE or joint_type == gs.JOINT_TYPE.PRISMATIC:
                     gear = -1  # TODO
-                    i_d = self.links_info[i_l].dof_start
+                    i_d = self.links_info[I_l].dof_start
                     self.dofs_state[i_d, i_b].act_length = gear * self.qpos[q_start, i_b]
                     self.dofs_state[i_d, i_b].qf_actuator = self.dofs_state[i_d, i_b].act_length
                 else:
-                    for i_d in range(self.links_info[i_l].dof_start, self.links_info[i_l].dof_end):
+                    for i_d in range(self.links_info[I_l].dof_start, self.links_info[I_l].dof_end):
                         self.dofs_state[i_d, i_b].act_length = 0.0
                         self.dofs_state[i_d, i_b].qf_actuator = self.dofs_state[i_d, i_b].act_length
 
@@ -2647,8 +2718,9 @@ def _func_bias_force(self):
             for i_b in range(self._B):
                 for i_l_ in range(self.n_awake_links[i_b]):
                     i_l = self.awake_links[i_l_, i_b]
+                    I_l = [i_l, i_b] if ti.static(self._options.batch_links_info) else i_l
 
-                    l_info = self.links_info[i_l]
+                    l_info = self.links_info[I_l]
                     dof_start = l_info.dof_start
                     dof_end = l_info.dof_end
 
@@ -2667,8 +2739,9 @@ def _func_bias_force(self):
         else:
             ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.ALL)
             for i_l, i_b in ti.ndrange(self.n_links, self._B):
+                I_l = [i_l, i_b] if ti.static(self._options.batch_links_info) else i_l
 
-                l_info = self.links_info[i_l]
+                l_info = self.links_info[I_l]
                 dof_start = l_info.dof_start
                 dof_end = l_info.dof_end
 
@@ -2727,10 +2800,11 @@ def _func_integrate(self):
             for i_b in range(self._B):
                 for i_l_ in range(self.n_awake_links[i_b]):
                     i_l = self.awake_links[i_l_, i_b]
-                    joint_type = self.links_info[i_l].joint_type
-                    dof_start = self.links_info[i_l].dof_start
-                    q_start = self.links_info[i_l].q_start
-                    q_end = self.links_info[i_l].q_end
+                    I_l = [i_l, i_b] if ti.static(self._options.batch_links_info) else i_l
+                    joint_type = self.links_info[I_l].joint_type
+                    dof_start = self.links_info[I_l].dof_start
+                    q_start = self.links_info[I_l].q_start
+                    q_end = self.links_info[I_l].q_end
 
                     if joint_type == gs.JOINT_TYPE.FREE:
                         rot = ti.Vector(
@@ -2795,10 +2869,11 @@ def _func_integrate(self):
 
             ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.ALL)
             for i_l, i_b in ti.ndrange(self.n_links, self._B):
-                joint_type = self.links_info[i_l].joint_type
-                dof_start = self.links_info[i_l].dof_start
-                q_start = self.links_info[i_l].q_start
-                q_end = self.links_info[i_l].q_end
+                I_l = [i_l, i_b] if ti.static(self._options.batch_links_info) else i_l
+                joint_type = self.links_info[I_l].joint_type
+                dof_start = self.links_info[I_l].dof_start
+                q_start = self.links_info[I_l].q_start
+                q_end = self.links_info[I_l].q_end
 
                 if joint_type == gs.JOINT_TYPE.FREE:
                     rot = ti.Vector(
@@ -2844,7 +2919,8 @@ def _func_integrate(self):
     def _func_integrate_dq_entity(self, dq, i_e, i_b, respect_joint_limit):
         e_info = self.entities_info[i_e]
         for i_l in range(e_info.link_start, e_info.link_end):
-            l_info = self.links_info[i_l]
+            I_l = [i_l, i_b] if ti.static(self._options.batch_links_info) else i_l
+            l_info = self.links_info[I_l]
             q_start = l_info.q_start
             dof_start = l_info.dof_start
             dq_start = l_info.dof_start - e_info.dof_start
@@ -2883,10 +2959,11 @@ def _func_integrate_dq_entity(self, dq, i_e, i_b, respect_joint_limit):
                     self.qpos[q_start + i_d_, i_b] = self.qpos[q_start + i_d_, i_b] + dq[dq_start + i_d_, i_b]
 
                     if respect_joint_limit:
+                        I_d = [dof_start + i_d_, i_b] if ti.static(self._options.batch_dofs_info) else dof_start + i_d_
                         self.qpos[q_start + i_d_, i_b] = ti.math.clamp(
                             self.qpos[q_start + i_d_, i_b],
-                            self.dofs_info[dof_start + i_d_].limit[0],
-                            self.dofs_info[dof_start + i_d_].limit[1],
+                            self.dofs_info[I_d].limit[0],
+                            self.dofs_info[I_d].limit[1],
                         )
 
     def substep_pre_coupling(self, f):
@@ -3206,12 +3283,13 @@ def _kernel_set_links_pos(
         ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.PARTIAL)
         for i_l_, i_b_ in ti.ndrange(links_idx.shape[0], envs_idx.shape[0]):
             i_l = links_idx[i_l_]
-            if self.links_info[i_l].is_fixed:  # change links_state directly as the link's pose is not contained in qpos
+            I_l = [i_l, i_b_] if ti.static(self._options.batch_links_info) else i_l
+            if self.links_info[I_l].is_fixed:  # change links_state directly as the link's pose is not contained in qpos
                 for i in ti.static(range(3)):
                     self.links_state[i_l, envs_idx[i_b_]].pos[i] = pos[i_b_, i_l_, i]
 
             else:  # free base link's pose is reflected in qpos, and links_state will be computed automatically
-                q_start = self.links_info[i_l].q_start
+                q_start = self.links_info[I_l].q_start
                 for i in ti.static(range(3)):
                     self.qpos[q_start + i, envs_idx[i_b_]] = pos[i_b_, i_l_, i]
 
@@ -3239,12 +3317,13 @@ def _kernel_set_links_quat(
         ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.PARTIAL)
         for i_l_, i_b_ in ti.ndrange(links_idx.shape[0], envs_idx.shape[0]):
             i_l = links_idx[i_l_]
-            if self.links_info[i_l].is_fixed:  # change links_state directly as the link's pose is not contained in qpos
+            I_l = [i_l, i_b_] if ti.static(self._options.batch_links_info) else i_l
+            if self.links_info[I_l].is_fixed:  # change links_state directly as the link's pose is not contained in qpos
                 for i in ti.static(range(4)):
                     self.links_state[i_l, envs_idx[i_b_]].quat[i] = quat[i_b_, i_l_, i]
 
             else:  # free base link's pose is reflected in qpos, and links_state will be computed automatically
-                q_start = self.links_info[i_l].q_start
+                q_start = self.links_info[I_l].q_start
                 for i in ti.static(range(4)):
                     self.qpos[q_start + i + 3, envs_idx[i_b_]] = quat[i_b_, i_l_, i]
 
@@ -3281,6 +3360,58 @@ def _kernel_set_links_COM_shift(
             for i in ti.static(range(3)):
                 self.links_state[links_idx[i_l_], envs_idx[i_b_]].i_pos_shift[i] = com[i_b_, i_l_, i]
 
+    def _set_links_info(self, tensor, links_idx, name, envs_idx=None):
+        if self._options.batch_links_info:
+            tensor, links_idx, envs_idx = self._validate_1D_io_variables(
+                tensor, links_idx, envs_idx, idx_name="links_idx"
+            )
+        else:
+            tensor, links_idx = self._validate_1D_io_variables(tensor, links_idx, idx_name="links_idx", batched=False)
+            envs_idx = torch.empty(())
+
+        if name == "invweight":
+            self._kernel_set_links_invweight(tensor, links_idx, envs_idx)
+        elif name == "inertial_mass":
+            self._kernel_set_links_inertial_mass(tensor, links_idx, envs_idx)
+        else:
+            gs.raise_exception(f"Invalid `name` {name}.")
+
+    def set_links_inertial_mass(self, invweight, links_idx, envs_idx=None):
+        self._set_links_info(invweight, links_idx, "inertial_mass", envs_idx)
+
+    @ti.kernel
+    def _kernel_set_links_inertial_mass(
+        self,
+        inertial_mass: ti.types.ndarray(),
+        links_idx: ti.types.ndarray(),
+        envs_idx: ti.types.ndarray(),
+    ):
+        ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.PARTIAL)
+        if ti.static(self._options.batch_links_info):
+            for i_l_, i_b_ in ti.ndrange(links_idx.shape[0], envs_idx.shape[0]):
+                self.links_info[links_idx[i_l_], envs_idx[i_b_]].inertial_mass = inertial_mass[i_b_, i_l_]
+        else:
+            for i_l_ in range(links_idx.shape[0]):
+                self.links_info[links_idx[i_l_]].inertial_mass = inertial_mass[i_l_]
+
+    def set_links_invweight(self, invweight, links_idx, envs_idx=None):
+        self._set_links_info(invweight, links_idx, "invweight", envs_idx)
+
+    @ti.kernel
+    def _kernel_set_links_invweight(
+        self,
+        invweight: ti.types.ndarray(),
+        links_idx: ti.types.ndarray(),
+        envs_idx: ti.types.ndarray(),
+    ):
+        ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.PARTIAL)
+        if ti.static(self._options.batch_links_info):
+            for i_l_, i_b_ in ti.ndrange(links_idx.shape[0], envs_idx.shape[0]):
+                self.links_info[links_idx[i_l_], envs_idx[i_b_]].invweight = invweight[i_b_, i_l_]
+        else:
+            for i_l_ in range(links_idx.shape[0]):
+                self.links_info[links_idx[i_l_]].invweight = invweight[i_l_]
+
     def set_geoms_friction_ratio(self, friction_ratio, geoms_idx, envs_idx=None):
         friction_ratio, geoms_idx, envs_idx = self._validate_1D_io_variables(
             friction_ratio, geoms_idx, envs_idx, idx_name="geoms_idx"
@@ -3337,59 +3468,195 @@ def _kernel_set_global_sol_params(self, sol_params: ti.types.ndarray()):
 
         ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.PARTIAL)
         for i, b in ti.ndrange(self.n_dofs, self._B):
+            I = [i, b] if ti.static(self._options.batch_dofs_info) else i
             for j in ti.static(range(7)):
-                self.dofs_info[i].sol_params[j] = sol_params[j]
+                self.dofs_info[I].sol_params[j] = sol_params[j]
 
-            self.dofs_info[i].sol_params[0] = self._substep_dt * 2
+            self.dofs_info[I].sol_params[0] = self._substep_dt * 2
 
-    def set_dofs_kp(self, kp, dofs_idx):
-        kp, dofs_idx = self._validate_1D_io_variables(kp, dofs_idx, batched=False)
-        self._kernel_set_dofs_kp(kp, dofs_idx)
+    def _set_dofs_info(self, tensor_list, dofs_idx, name, envs_idx=None):
+        if self._options.batch_dofs_info:
+            for i, tensor in enumerate(tensor_list):
+                if i == (len(tensor_list) - 1):
+                    tensor_list[i], dofs_idx, envs_idx = self._validate_1D_io_variables(tensor, dofs_idx, envs_idx)
+                else:
+                    tensor_list[i], _, _ = self._validate_1D_io_variables(tensor, dofs_idx, envs_idx)
+        else:
+            for i, tensor in enumerate(tensor_list):
+                if i == (len(tensor_list) - 1):
+                    tensor_list[i], _ = self._validate_1D_io_variables(tensor, dofs_idx, batched=False)
+                else:
+                    tensor_list[i], dofs_idx = self._validate_1D_io_variables(tensor, dofs_idx, batched=False)
+            envs_idx = torch.empty(())
+
+        if name == "kp":
+            self._kernel_set_dofs_kp(tensor_list[0], dofs_idx, envs_idx)
+        elif name == "kv":
+            self._kernel_set_dofs_kv(tensor_list[0], dofs_idx, envs_idx)
+        elif name == "force_range":
+            self._kernel_set_dofs_force_range(tensor_list[0], tensor_list[1], dofs_idx, envs_idx)
+        elif name == "stiffness":
+            self._kernel_set_dofs_stiffness(tensor_list[0], dofs_idx, envs_idx)
+        elif name == "invweight":
+            self._kernel_set_dofs_invweight(tensor_list[0], dofs_idx, envs_idx)
+        elif name == "armature":
+            self._kernel_set_dofs_armature(tensor_list[0], dofs_idx, envs_idx)
+        elif name == "damping":
+            self._kernel_set_dofs_damping(tensor_list[0], dofs_idx, envs_idx)
+        elif name == "limit":
+            self._kernel_set_dofs_limit(tensor_list[0], tensor_list[1], dofs_idx, envs_idx)
+        else:
+            gs.raise_exception(f"Invalid `name` {name}.")
+
+    def set_dofs_kp(self, kp, dofs_idx, envs_idx=None):
+        self._set_dofs_info([kp], dofs_idx, "kp", envs_idx)
 
     @ti.kernel
     def _kernel_set_dofs_kp(
         self,
         kp: ti.types.ndarray(),
         dofs_idx: ti.types.ndarray(),
+        envs_idx: ti.types.ndarray(),
     ):
         ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.PARTIAL)
-        for i_d_ in range(dofs_idx.shape[0]):
-            self.dofs_info[dofs_idx[i_d_]].kp = kp[i_d_]
+        if ti.static(self._options.batch_dofs_info):
+            for i_d_, i_b_ in ti.ndrange(dofs_idx.shape[0], envs_idx.shape[0]):
+                self.dofs_info[dofs_idx[i_d_], envs_idx[i_b_]].kp = kp[i_b_, i_d_]
+        else:
+            for i_d_ in range(dofs_idx.shape[0]):
+                self.dofs_info[dofs_idx[i_d_]].kp = kp[i_d_]
 
-    def set_dofs_kv(self, kv, dofs_idx):
-        kv, dofs_idx = self._validate_1D_io_variables(kv, dofs_idx, batched=False)
-        self._kernel_set_dofs_kv(kv, dofs_idx)
+    def set_dofs_kv(self, kv, dofs_idx, envs_idx=None):
+        self._set_dofs_info([kv], dofs_idx, "kv", envs_idx)
 
     @ti.kernel
     def _kernel_set_dofs_kv(
         self,
         kv: ti.types.ndarray(),
         dofs_idx: ti.types.ndarray(),
+        envs_idx: ti.types.ndarray(),
     ):
         ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.PARTIAL)
-        for i_d_ in range(dofs_idx.shape[0]):
-            self.dofs_info[dofs_idx[i_d_]].kv = kv[i_d_]
+        if ti.static(self._options.batch_dofs_info):
+            for i_d_, i_b_ in ti.ndrange(dofs_idx.shape[0], envs_idx.shape[0]):
+                self.dofs_info[dofs_idx[i_d_], envs_idx[i_b_]].kv = kv[i_b_, i_d_]
+        else:
+            for i_d_ in range(dofs_idx.shape[0]):
+                self.dofs_info[dofs_idx[i_d_]].kv = kv[i_d_]
 
-    def set_dofs_force_range(self, lower, upper, dofs_idx):
-        lower, _ = self._validate_1D_io_variables(lower, dofs_idx, batched=False)
-        upper, dofs_idx = self._validate_1D_io_variables(upper, dofs_idx, batched=False)
+    def set_dofs_force_range(self, lower, upper, dofs_idx, envs_idx=None):
+        self._set_dofs_info([lower, upper], dofs_idx, "force_range", envs_idx)
 
-        if (lower > upper).any():
-            gs.raise_exception("`lower` should be less than or equal to `upper`.")
+    @ti.kernel
+    def _kernel_set_dofs_force_range(
+        self,
+        lower: ti.types.ndarray(),
+        upper: ti.types.ndarray(),
+        dofs_idx: ti.types.ndarray(),
+        envs_idx: ti.types.ndarray(),
+    ):
+        ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.PARTIAL)
+        if ti.static(self._options.batch_dofs_info):
+            for i_d_, i_b_ in ti.ndrange(dofs_idx.shape[0], envs_idx.shape[0]):
+                self.dofs_info[dofs_idx[i_d_], envs_idx[i_b_]].force_range[0] = lower[i_b_, i_d_]
+                self.dofs_info[dofs_idx[i_d_], envs_idx[i_b_]].force_range[1] = upper[i_b_, i_d_]
+        else:
+            for i_d_ in range(dofs_idx.shape[0]):
+                self.dofs_info[dofs_idx[i_d_]].force_range[0] = lower[i_d_]
+                self.dofs_info[dofs_idx[i_d_]].force_range[1] = upper[i_d_]
 
-        self._kernel_set_dofs_force_range(lower, upper, dofs_idx)
+    def set_dofs_stiffness(self, stiffness, dofs_idx, envs_idx=None):
+        self._set_dofs_info([stiffness], dofs_idx, "stiffness", envs_idx)
 
     @ti.kernel
-    def _kernel_set_dofs_force_range(
+    def _kernel_set_dofs_stiffness(
+        self,
+        stiffness: ti.types.ndarray(),
+        dofs_idx: ti.types.ndarray(),
+        envs_idx: ti.types.ndarray(),
+    ):
+        ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.PARTIAL)
+        if ti.static(self._options.batch_dofs_info):
+            for i_d_, i_b_ in ti.ndrange(dofs_idx.shape[0], envs_idx.shape[0]):
+                self.dofs_info[dofs_idx[i_d_], envs_idx[i_b_]].stiffness = stiffness[i_b_, i_d_]
+        else:
+            for i_d_ in range(dofs_idx.shape[0]):
+                self.dofs_info[dofs_idx[i_d_]].stiffness = stiffness[i_d_]
+
+    def set_dofs_invweight(self, invweight, dofs_idx, envs_idx=None):
+        self._set_dofs_info([invweight], dofs_idx, "invweight", envs_idx)
+
+    @ti.kernel
+    def _kernel_set_dofs_invweight(
+        self,
+        invweight: ti.types.ndarray(),
+        dofs_idx: ti.types.ndarray(),
+        envs_idx: ti.types.ndarray(),
+    ):
+        ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.PARTIAL)
+        if ti.static(self._options.batch_dofs_info):
+            for i_d_, i_b_ in ti.ndrange(dofs_idx.shape[0], envs_idx.shape[0]):
+                self.dofs_info[dofs_idx[i_d_], envs_idx[i_b_]].invweight = invweight[i_b_, i_d_]
+        else:
+            for i_d_ in range(dofs_idx.shape[0]):
+                self.dofs_info[dofs_idx[i_d_]].invweight = invweight[i_d_]
+
+    def set_dofs_armature(self, armature, dofs_idx, envs_idx=None):
+        self._set_dofs_info([armature], dofs_idx, "armature", envs_idx)
+
+    @ti.kernel
+    def _kernel_set_dofs_armature(
+        self,
+        armature: ti.types.ndarray(),
+        dofs_idx: ti.types.ndarray(),
+        envs_idx: ti.types.ndarray(),
+    ):
+        ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.PARTIAL)
+        if ti.static(self._options.batch_dofs_info):
+            for i_d_, i_b_ in ti.ndrange(dofs_idx.shape[0], envs_idx.shape[0]):
+                self.dofs_info[dofs_idx[i_d_], envs_idx[i_b_]].armature = armature[i_b_, i_d_]
+        else:
+            for i_d_ in range(dofs_idx.shape[0]):
+                self.dofs_info[dofs_idx[i_d_]].armature = armature[i_d_]
+
+    def set_dofs_damping(self, damping, dofs_idx, envs_idx=None):
+        self._set_dofs_info([damping], dofs_idx, "damping", envs_idx)
+
+    @ti.kernel
+    def _kernel_set_dofs_damping(
+        self,
+        damping: ti.types.ndarray(),
+        dofs_idx: ti.types.ndarray(),
+        envs_idx: ti.types.ndarray(),
+    ):
+        ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.PARTIAL)
+        if ti.static(self._options.batch_dofs_info):
+            for i_d_, i_b_ in ti.ndrange(dofs_idx.shape[0], envs_idx.shape[0]):
+                self.dofs_info[dofs_idx[i_d_], envs_idx[i_b_]].damping = damping[i_b_, i_d_]
+        else:
+            for i_d_ in range(dofs_idx.shape[0]):
+                self.dofs_info[dofs_idx[i_d_]].damping = damping[i_d_]
+
+    def set_dofs_limit(self, lower, upper, dofs_idx, envs_idx=None):
+        self._set_dofs_info([lower, upper], dofs_idx, "limit", envs_idx)
+
+    @ti.kernel
+    def _kernel_set_dofs_limit(
         self,
         lower: ti.types.ndarray(),
         upper: ti.types.ndarray(),
         dofs_idx: ti.types.ndarray(),
+        envs_idx: ti.types.ndarray(),
     ):
         ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.PARTIAL)
-        for i_d_ in range(dofs_idx.shape[0]):
-            self.dofs_info[dofs_idx[i_d_]].force_range[0] = lower[i_d_]
-            self.dofs_info[dofs_idx[i_d_]].force_range[1] = upper[i_d_]
+        if ti.static(self._options.batch_dofs_info):
+            for i_d_, i_b_ in ti.ndrange(dofs_idx.shape[0], envs_idx.shape[0]):
+                self.dofs_info[dofs_idx[i_d_], envs_idx[i_b_]].limit[0] = lower[i_b_, i_d_]
+                self.dofs_info[dofs_idx[i_d_], envs_idx[i_b_]].limit[1] = upper[i_b_, i_d_]
+        else:
+            for i_d_ in range(dofs_idx.shape[0]):
+                self.dofs_info[dofs_idx[i_d_]].limit[0] = lower[i_d_]
+                self.dofs_info[dofs_idx[i_d_]].limit[1] = upper[i_d_]
 
     def set_dofs_velocity(self, velocity, dofs_idx, envs_idx=None):
         velocity, dofs_idx, envs_idx = self._validate_1D_io_variables(velocity, dofs_idx, envs_idx)
@@ -3437,7 +3704,8 @@ def _kernel_set_dofs_position(
         for i_e, i_b_ in ti.ndrange(self.n_entities, envs_idx.shape[0]):
             i_b = envs_idx[i_b_]
             for i_l in range(self.entities_info[i_e].link_start, self.entities_info[i_e].link_end):
-                l_info = self.links_info[i_l]
+                I_l = [i_l, i_b] if ti.static(self._options.batch_links_info) else i_l
+                l_info = self.links_info[I_l]
 
                 if l_info.joint_type == gs.JOINT_TYPE.FREE:
                     xyz = ti.Vector(
@@ -3664,6 +3932,60 @@ def _kernel_get_links_COM_shift(
             for i in ti.static(range(3)):
                 tensor[i_b_, i_l_, i] = self.links_state[links_idx[i_l_], envs_idx[i_b_]].i_pos_shift[i]
 
+    def _get_links_info(self, links_idx, name, envs_idx=None):
+        if self._options.batch_links_info:
+            tensor, links_idx, envs_idx = self._validate_1D_io_variables(
+                None, links_idx, envs_idx, idx_name="links_idx"
+            )
+        else:
+            tensor, links_idx = self._validate_1D_io_variables(None, links_idx, idx_name="links_idx", batched=False)
+            envs_idx = torch.empty(())
+
+        if name == "invweight":
+            self._kernel_get_links_invweight(tensor, links_idx, envs_idx)
+            return tensor
+        elif name == "inertial_mass":
+            self._kernel_get_links_inertial_mass(tensor, links_idx, envs_idx)
+            return tensor
+        else:
+            gs.raise_exception(f"Invalid `name` {name}.")
+
+    def get_links_inertial_mass(self, links_idx, envs_idx=None):
+        return self._get_links_info(links_idx, "inertial_mass", envs_idx)
+
+    @ti.kernel
+    def _kernel_get_links_inertial_mass(
+        self,
+        tensor: ti.types.ndarray(),
+        links_idx: ti.types.ndarray(),
+        envs_idx: ti.types.ndarray(),
+    ):
+        ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.PARTIAL)
+        if ti.static(self._options.batch_links_info):
+            for i_l_, i_b_ in ti.ndrange(links_idx.shape[0], envs_idx.shape[0]):
+                tensor[i_b_, i_l_] = self.links_info[links_idx[i_l_], envs_idx[i_b_]].inertial_mass
+        else:
+            for i_l_ in range(links_idx.shape[0]):
+                tensor[i_l_] = self.links_info[links_idx[i_l_]].inertial_mass
+
+    def get_links_invweight(self, links_idx, envs_idx=None):
+        return self._get_links_info(links_idx, "invweight", envs_idx)
+
+    @ti.kernel
+    def _kernel_get_links_invweight(
+        self,
+        tensor: ti.types.ndarray(),
+        links_idx: ti.types.ndarray(),
+        envs_idx: ti.types.ndarray(),
+    ):
+        ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.PARTIAL)
+        if ti.static(self._options.batch_links_info):
+            for i_l_, i_b_ in ti.ndrange(links_idx.shape[0], envs_idx.shape[0]):
+                tensor[i_b_, i_l_] = self.links_info[links_idx[i_l_], envs_idx[i_b_]].invweight
+        else:
+            for i_l_ in range(links_idx.shape[0]):
+                tensor[i_l_] = self.links_info[links_idx[i_l_]].invweight
+
     def get_geoms_friction_ratio(self, geoms_idx, envs_idx=None):
         tensor, geoms_idx, envs_idx = self._validate_1D_io_variables(None, geoms_idx, envs_idx, idx_name="geoms_idx")
 
@@ -3767,20 +4089,21 @@ def _kernel_get_dofs_control_force(
         for i_d_, i_b_ in ti.ndrange(dofs_idx.shape[0], envs_idx.shape[0]):
             i_d = dofs_idx[i_d_]
             i_b = envs_idx[i_b_]
+            I_d = [i_d, i_b] if ti.static(self._options.batch_dofs_info) else i_d
             force = gs.ti_float(0.0)
             if self.dofs_state[i_d, i_b].ctrl_mode == gs.CTRL_MODE.FORCE:
                 force = self.dofs_state[i_d, i_b].ctrl_force
             elif self.dofs_state[i_d, i_b].ctrl_mode == gs.CTRL_MODE.VELOCITY:
-                force = self.dofs_info[i_d].kv * (self.dofs_state[i_d, i_b].ctrl_vel - self.dofs_state[i_d, i_b].vel)
+                force = self.dofs_info[I_d].kv * (self.dofs_state[i_d, i_b].ctrl_vel - self.dofs_state[i_d, i_b].vel)
             elif self.dofs_state[i_d, i_b].ctrl_mode == gs.CTRL_MODE.POSITION:
                 force = (
-                    self.dofs_info[i_d].kp * (self.dofs_state[i_d, i_b].ctrl_pos - self.dofs_state[i_d, i_b].pos)
-                    - self.dofs_info[i_d].kv * self.dofs_state[i_d, i_b].vel
+                    self.dofs_info[I_d].kp * (self.dofs_state[i_d, i_b].ctrl_pos - self.dofs_state[i_d, i_b].pos)
+                    - self.dofs_info[I_d].kv * self.dofs_state[i_d, i_b].vel
                 )
             tensor[i_b_, i_d_] = ti.math.clamp(
                 force,
-                self.dofs_info[i_d].force_range[0],
-                self.dofs_info[i_d].force_range[1],
+                self.dofs_info[I_d].force_range[0],
+                self.dofs_info[I_d].force_range[1],
             )
 
     @ti.kernel
@@ -3828,29 +4151,61 @@ def get_dofs_force_range(self, dofs_idx, envs_idx=None):
     def get_dofs_limit(self, dofs_idx, envs_idx=None):
         return self._get_dofs_info(dofs_idx, "limit", envs_idx)
 
+    def get_dofs_stiffness(self, dofs_idx, envs_idx=None):
+        return self._get_dofs_info(dofs_idx, "stiffness", envs_idx)
+
+    def get_dofs_invweight(self, dofs_idx, envs_idx=None):
+        return self._get_dofs_info(dofs_idx, "invweight", envs_idx)
+
+    def get_dofs_armature(self, dofs_idx, envs_idx=None):
+        return self._get_dofs_info(dofs_idx, "armature", envs_idx)
+
+    def get_dofs_damping(self, dofs_idx, envs_idx=None):
+        return self._get_dofs_info(dofs_idx, "damping", envs_idx)
+
     def _get_dofs_info(self, dofs_idx, name, envs_idx=None):
-        tensor, dofs_idx = self._validate_1D_io_variables(None, dofs_idx, batched=False)
+        if self._options.batch_dofs_info:
+            tensor, dofs_idx, envs_idx = self._validate_1D_io_variables(None, dofs_idx, envs_idx)
+        else:
+            tensor, dofs_idx = self._validate_1D_io_variables(None, dofs_idx, batched=False)
+            envs_idx = torch.empty(())
 
         if name == "kp":
-            self._kernel_get_dofs_kp(tensor, dofs_idx)
+            self._kernel_get_dofs_kp(tensor, dofs_idx, envs_idx)
             return tensor
 
         elif name == "kv":
-            self._kernel_get_dofs_kv(tensor, dofs_idx)
+            self._kernel_get_dofs_kv(tensor, dofs_idx, envs_idx)
             return tensor
 
         elif name == "force_range":
             lower = torch.empty_like(tensor)
             upper = torch.empty_like(tensor)
-            self._kernel_get_dofs_force_range(lower, upper, dofs_idx)
+            self._kernel_get_dofs_force_range(lower, upper, dofs_idx, envs_idx)
             return lower, upper
 
         elif name == "limit":
             lower = torch.empty_like(tensor)
             upper = torch.empty_like(tensor)
-            self._kernel_get_dofs_limit(lower, upper, dofs_idx)
+            self._kernel_get_dofs_limit(lower, upper, dofs_idx, envs_idx)
             return lower, upper
 
+        elif name == "stiffness":
+            self._kernel_get_dofs_stiffness(tensor, dofs_idx, envs_idx)
+            return tensor
+
+        elif name == "invweight":
+            self._kernel_get_dofs_invweight(tensor, dofs_idx, envs_idx)
+            return tensor
+
+        elif name == "armature":
+            self._kernel_get_dofs_armature(tensor, dofs_idx, envs_idx)
+            return tensor
+
+        elif name == "damping":
+            self._kernel_get_dofs_damping(tensor, dofs_idx, envs_idx)
+            return tensor
+
         else:
             gs.raise_exception()
 
@@ -3859,20 +4214,30 @@ def _kernel_get_dofs_kp(
         self,
         tensor: ti.types.ndarray(),
         dofs_idx: ti.types.ndarray(),
+        envs_idx: ti.types.ndarray(),
     ):
         ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.PARTIAL)
-        for i_d_ in range(dofs_idx.shape[0]):
-            tensor[i_d_] = self.dofs_info[dofs_idx[i_d_]].kp
+        if ti.static(self._options.batch_dofs_info):
+            for i_d_, i_b_ in ti.ndrange(dofs_idx.shape[0], envs_idx.shape[0]):
+                tensor[i_b_, i_d_] = self.dofs_info[dofs_idx[i_d_], envs_idx[i_b_]].kp
+        else:
+            for i_d_ in range(dofs_idx.shape[0]):
+                tensor[i_d_] = self.dofs_info[dofs_idx[i_d_]].kp
 
     @ti.kernel
     def _kernel_get_dofs_kv(
         self,
         tensor: ti.types.ndarray(),
         dofs_idx: ti.types.ndarray(),
+        envs_idx: ti.types.ndarray(),
     ):
         ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.PARTIAL)
-        for i_d_ in range(dofs_idx.shape[0]):
-            tensor[i_d_] = self.dofs_info[dofs_idx[i_d_]].kv
+        if ti.static(self._options.batch_dofs_info):
+            for i_d_, i_b_ in ti.ndrange(dofs_idx.shape[0], envs_idx.shape[0]):
+                tensor[i_b_, i_d_] = self.dofs_info[dofs_idx[i_d_], envs_idx[i_b_]].kv
+        else:
+            for i_d_ in range(dofs_idx.shape[0]):
+                tensor[i_d_] = self.dofs_info[dofs_idx[i_d_]].kv
 
     @ti.kernel
     def _kernel_get_dofs_force_range(
@@ -3880,11 +4245,17 @@ def _kernel_get_dofs_force_range(
         lower: ti.types.ndarray(),
         upper: ti.types.ndarray(),
         dofs_idx: ti.types.ndarray(),
+        envs_idx: ti.types.ndarray(),
     ):
         ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.PARTIAL)
-        for i_d_ in range(dofs_idx.shape[0]):
-            lower[i_d_] = self.dofs_info[dofs_idx[i_d_]].force_range[0]
-            upper[i_d_] = self.dofs_info[dofs_idx[i_d_]].force_range[1]
+        if ti.static(self._options.batch_dofs_info):
+            for i_d_, i_b_ in ti.ndrange(dofs_idx.shape[0], envs_idx.shape[0]):
+                lower[i_b_, i_d_] = self.dofs_info[dofs_idx[i_d_], envs_idx[i_b_]].force_range[0]
+                upper[i_b_, i_d_] = self.dofs_info[dofs_idx[i_d_], envs_idx[i_b_]].force_range[1]
+        else:
+            for i_d_ in range(dofs_idx.shape[0]):
+                lower[i_d_] = self.dofs_info[dofs_idx[i_d_]].force_range[0]
+                upper[i_d_] = self.dofs_info[dofs_idx[i_d_]].force_range[1]
 
     @ti.kernel
     def _kernel_get_dofs_limit(
@@ -3892,11 +4263,77 @@ def _kernel_get_dofs_limit(
         lower: ti.types.ndarray(),
         upper: ti.types.ndarray(),
         dofs_idx: ti.types.ndarray(),
+        envs_idx: ti.types.ndarray(),
     ):
         ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.PARTIAL)
-        for i_d_ in range(dofs_idx.shape[0]):
-            lower[i_d_] = self.dofs_info[dofs_idx[i_d_]].limit[0]
-            upper[i_d_] = self.dofs_info[dofs_idx[i_d_]].limit[1]
+        if ti.static(self._options.batch_dofs_info):
+            for i_d_, i_b_ in ti.ndrange(dofs_idx.shape[0], envs_idx.shape[0]):
+                lower[i_b_, i_d_] = self.dofs_info[dofs_idx[i_d_], envs_idx[i_b_]].limit[0]
+                upper[i_b_, i_d_] = self.dofs_info[dofs_idx[i_d_], envs_idx[i_b_]].limit[1]
+        else:
+            for i_d_ in range(dofs_idx.shape[0]):
+                lower[i_d_] = self.dofs_info[dofs_idx[i_d_]].limit[0]
+                upper[i_d_] = self.dofs_info[dofs_idx[i_d_]].limit[1]
+
+    @ti.kernel
+    def _kernel_get_dofs_stiffness(
+        self,
+        tensor: ti.types.ndarray(),
+        dofs_idx: ti.types.ndarray(),
+        envs_idx: ti.types.ndarray(),
+    ):
+        ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.PARTIAL)
+        if ti.static(self._options.batch_dofs_info):
+            for i_d_, i_b_ in ti.ndrange(dofs_idx.shape[0], envs_idx.shape[0]):
+                tensor[i_b_, i_d_] = self.dofs_info[dofs_idx[i_d_], envs_idx[i_b_]].stiffness
+        else:
+            for i_d_ in range(dofs_idx.shape[0]):
+                tensor[i_d_] = self.dofs_info[dofs_idx[i_d_]].stiffness
+
+    @ti.kernel
+    def _kernel_get_dofs_invweight(
+        self,
+        tensor: ti.types.ndarray(),
+        dofs_idx: ti.types.ndarray(),
+        envs_idx: ti.types.ndarray(),
+    ):
+        ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.PARTIAL)
+        if ti.static(self._options.batch_dofs_info):
+            for i_d_, i_b_ in ti.ndrange(dofs_idx.shape[0], envs_idx.shape[0]):
+                tensor[i_b_, i_d_] = self.dofs_info[dofs_idx[i_d_], envs_idx[i_b_]].invweight
+        else:
+            for i_d_ in range(dofs_idx.shape[0]):
+                tensor[i_d_] = self.dofs_info[dofs_idx[i_d_]].invweight
+
+    @ti.kernel
+    def _kernel_get_dofs_armature(
+        self,
+        tensor: ti.types.ndarray(),
+        dofs_idx: ti.types.ndarray(),
+        envs_idx: ti.types.ndarray(),
+    ):
+        ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.PARTIAL)
+        if ti.static(self._options.batch_dofs_info):
+            for i_d_, i_b_ in ti.ndrange(dofs_idx.shape[0], envs_idx.shape[0]):
+                tensor[i_b_, i_d_] = self.dofs_info[dofs_idx[i_d_], envs_idx[i_b_]].armature
+        else:
+            for i_d_ in range(dofs_idx.shape[0]):
+                tensor[i_d_] = self.dofs_info[dofs_idx[i_d_]].armature
+
+    @ti.kernel
+    def _kernel_get_dofs_damping(
+        self,
+        tensor: ti.types.ndarray(),
+        dofs_idx: ti.types.ndarray(),
+        envs_idx: ti.types.ndarray(),
+    ):
+        ti.loop_config(serialize=self._para_level < gs.PARA_LEVEL.PARTIAL)
+        if ti.static(self._options.batch_dofs_info):
+            for i_d_, i_b_ in ti.ndrange(dofs_idx.shape[0], envs_idx.shape[0]):
+                tensor[i_b_, i_d_] = self.dofs_info[dofs_idx[i_d_], envs_idx[i_b_]].damping
+        else:
+            for i_d_ in range(dofs_idx.shape[0]):
+                tensor[i_d_] = self.dofs_info[dofs_idx[i_d_]].damping
 
     @ti.kernel
     def _kernel_set_drone_rpm(
diff --git a/genesis/ext/pyrender/viewer.py b/genesis/ext/pyrender/viewer.py
index 7b6c8700..e2c20ad9 100644
--- a/genesis/ext/pyrender/viewer.py
+++ b/genesis/ext/pyrender/viewer.py
@@ -14,15 +14,14 @@
 
 import genesis as gs
 
+from tkinter import Tk
+from tkinter import filedialog
+
 try:
-    from Tkinter import Tk
-    from Tkinter import tkFileDialog as filedialog
-except Exception:
-    try:
-        from tkinter import Tk
-        from tkinter import filedialog as filedialog
-    except Exception:
-        pass
+    root = Tk()
+    root.withdraw()
+except:
+    pass
 
 import pyglet
 from moviepy.video.io.ffmpeg_writer import FFMPEG_VideoWriter
@@ -985,7 +984,6 @@ def _get_save_filename(self, file_exts):
         }
         filetypes = [file_types[x] for x in file_exts]
         try:
-            root = Tk()
             save_dir = self.viewer_flags["save_directory"]
             if save_dir is None:
                 save_dir = os.getcwd()
@@ -995,7 +993,6 @@ def _get_save_filename(self, file_exts):
         except Exception:
             return None
 
-        root.destroy()
         if filename == ():
             return None
         return filename
diff --git a/genesis/options/solvers.py b/genesis/options/solvers.py
index 54e80b18..1b3b6c51 100644
--- a/genesis/options/solvers.py
+++ b/genesis/options/solvers.py
@@ -185,6 +185,10 @@ class RigidOptions(Options):
     integrator: gs.integrator = gs.integrator.approximate_implicitfast
     IK_max_targets: int = 6
 
+    # batching info
+    batch_links_info: Optional[bool] = False
+    batch_dofs_info: Optional[bool] = False
+
     # constraint solver
     constraint_solver: gs.constraint_solver = gs.constraint_solver.CG
     iterations: int = 100
diff --git a/genesis/repr_base.py b/genesis/repr_base.py
index 86f18cf0..b81d0e07 100644
--- a/genesis/repr_base.py
+++ b/genesis/repr_base.py
@@ -1,3 +1,5 @@
+import inspect
+
 import genesis as gs
 import genesis.utils.repr as ru
 from genesis.styles import colors, formats, styless
@@ -31,9 +33,16 @@ def _repr_brief(self):
             repr_str += f", material: {self.material}"
         return repr_str
 
+    def _is_debugger(self) -> bool:
+        """Detect if running under a debugger (VSCode or PyCharm)."""
+        for frame in inspect.stack():
+            if any(module in frame.filename for module in ("debugpy", "ptvsd", "pydevd")):
+                return True
+        return False
+
     def __repr__(self):
-        if not __debug__:
-            self.__colorized__repr__()
+        if not self._is_debugger():
+            return self.__colorized__repr__()
 
     def __colorized__repr__(self) -> str:
         all_attrs = self.__dir__()
diff --git a/genesis/utils/mjcf.py b/genesis/utils/mjcf.py
index b961f288..45528816 100644
--- a/genesis/utils/mjcf.py
+++ b/genesis/utils/mjcf.py
@@ -45,7 +45,6 @@ def parse_mjcf(path):
 
 
 def parse_link(mj, i_l, q_offset, dof_offset, scale):
-
     # mj.body
     l_info = dict()
 
@@ -286,6 +285,32 @@ def parse_geom(mj, i_g, scale, convexify, surface, xml_path):
         tmesh = trimesh.creation.box(extents=mj.geom_size[i_g, :3] * 2)
         gs_type = gs.GEOM_TYPE.BOX
 
+        # TODO: not sure if it is the right way to load texture for box
+        mat_id = mj.geom_matid[i_g]
+        if mat_id >= 0:
+            mat_id = mj.geom_matid[i_g]
+            tex_id = next((x for x in mj.mat_texid[mat_id] if x != -1), None)
+
+            if tex_id is not None:
+                tex_path = mj.paths[mj.tex_pathadr[tex_id] :].decode("utf-8").split("\x00")[0]
+                texturedir = extract_compiler_attributes(xml_path)["texturedir"]
+                assets_dir = os.path.join(get_assets_dir(), os.path.join(os.path.dirname(xml_path), texturedir))
+
+                uv_coordinates = tmesh.vertices[:, :2].copy()
+                uv_coordinates -= uv_coordinates.min(axis=0)
+                uv_coordinates /= uv_coordinates.max(axis=0)
+                image = Image.open(os.path.join(assets_dir, tex_path))
+                image_array = np.array(image)
+                if image_array.ndim == 2:  # convert gray image to RGBA
+                    rgba = np.zeros((image_array.shape[0], image_array.shape[1], 4), dtype=np.uint8)
+                    rgba[:, :, :3] = image_array[:, :, None]
+                    rgba[:, :, 3] = 255
+                    image_array = rgba
+                tex_repeat = mj.mat_texrepeat[mat_id].astype(int)
+                image_array = np.tile(image_array, (tex_repeat[0], tex_repeat[1], 1))
+                visual = TextureVisuals(uv=uv_coordinates, image=Image.fromarray(image_array))
+                tmesh.visual = visual
+
     elif mj_type == mujoco.mjtGeom.mjGEOM_MESH:
         i = mj.geom_dataid[i_g]
         last = (i + 1) >= mj.nmesh