Make Rapier optional and support Rapier f64

CHANGELOG.md

@@ -1,4 +1,8 @@
 # Unreleased
+- Make Rapier an optional dependency. Enable either the `rapier-f32` or `rapier-f64`
+  feature for coupling with Rapier rigid bodies and colliders support.
+- Add support for the `f64` version of Rapier with the `rapier-f64` feature.
+  The MPM simulation still uses single-precision values.
 - Update to Rapier 0.32. This migrates most public APIs and internals to use `glam` instead of `nalgebra`.
 - Fix a GPU validation error / panic on simulations with more than ~4.19M particles, caused by
   compute kernels dispatching more than 65535 workgroups along a single dimension. The affected 

crates/slosh2d/Cargo.toml

@@ -16,11 +16,15 @@ path = "../../src/lib.rs"
 required-features = ["dim2"]
 
 [features]
-default = ["dim2", "cpic", "node_particle_lists"]
+default = ["dim2", "rapier-f32", "cpic", "node_particle_lists"]
 dim2 = []
 cpic = []
 node_particle_lists = []
 
+rapier-f32 = ["dep:rapier2d", "rapier"]
+rapier-f64 = ["dep:rapier2d-f64", "rapier"]
+rapier = []
+
 comptime = ["slosh_testbed2d/comptime", "stensor/comptime"]
 runtime = ["slosh_testbed2d/runtime", "stensor/runtime"]
 webgpu = ["slosh_testbed2d/webgpu", "stensor/webgpu"]
@@ -38,8 +42,8 @@ static_assertions = { workspace = true }
 #wgebra = "0.2"
 bvh = "0.12.0"
 
-# TODO: make rapier optional?
-rapier2d = "0.32"
+rapier2d = { version = "0.32", optional = true }
+rapier2d-f64 = { version = "0.32", optional = true }
 
 # For wasm?
 getrandom = { version = "0.3.1", features = ["wasm_js"] }

crates/slosh3d/Cargo.toml

@@ -16,11 +16,15 @@ path = "../../src/lib.rs"
 required-features = ["dim3"]
 
 [features]
-default = ["dim3", "cpic", "node_particle_lists"]
+default = ["dim3", "rapier-f32", "cpic", "node_particle_lists"]
 dim3 = []
 cpic = []
 node_particle_lists = []
 
+rapier-f32 = ["dep:rapier3d", "rapier"]
+rapier-f64 = ["dep:rapier3d-f64", "rapier"]
+rapier = []
+
 comptime = ["stensor/comptime", "slosh_testbed3d/comptime"]
 runtime = ["stensor/runtime", "slosh_testbed3d/runtime"]
 webgpu = ["stensor/webgpu", "slosh_testbed3d/webgpu"]
@@ -38,8 +42,8 @@ bvh = "0.12.0"
 
 serde = "1"
 
-# TODO: make rapier optional?
-rapier3d = "0.32"
+rapier3d = { version = "0.32", optional = true }
+rapier3d-f64 = { version = "0.32", optional = true }
 
 [dev-dependencies]
 futures-test = "0.3"

crates/slosh_testbed2d/Cargo.toml

@@ -38,7 +38,7 @@ futures-test = "0.3"
 serial_test = "3"
 approx = "0.5"
 async-std = { version = "1", features = ["attributes"] }
-slosh2d = { version = "0.7", path = "../slosh2d" }
+slosh2d = { version = "0.7", path = "../slosh2d", features = ["rapier-f32"] }
 regex = "1"
 web-time = "1"
 

crates/slosh_testbed3d/Cargo.toml

@@ -38,7 +38,7 @@ futures-test = "0.3"
 serial_test = "3"
 approx = "0.5"
 async-std = { version = "1", features = ["attributes"] }
-slosh3d = { version = "0.7", path = "../slosh3d" }
+slosh3d = { version = "0.7", path = "../slosh3d", features = ["rapier-f32"] }
 regex = "1"
 web-time = "1"
 

src/lib.rs

@@ -18,7 +18,9 @@
 //! # Features
 //!
 //! - `dim2`: Enable 2D simulation mode (mutually exclusive with `dim3`)
-//! - `dim3`: Enable 3D simulation mode (default, mutually exclusive with `dim2`)
+//! - `dim3`: Enable 3D simulation mode (mutually exclusive with `dim2`)
+//! - `rapier-f32`: Support coupling with Rapier using `f32` precision (mutually exclusive with `rapier-f64`)
+//! - `rapier-f64`: Support coupling with Rapier using `f64` precision (mutually exclusive with `rapier-f32`)
 //!
 //! # Example
 //!
@@ -55,10 +57,24 @@
 #![allow(clippy::module_inception)]
 #![allow(missing_docs)]
 
-#[cfg(feature = "dim2")]
+#[cfg(all(feature = "rapier-f32", feature = "rapier-f64"))]
+compile_error!(
+    "Features `rapier-f32` and `rapier-f64` are mutually exclusive. Please enable only one of them."
+);
+#[cfg(all(
+    feature = "rapier",
+    not(any(feature = "rapier-f32", feature = "rapier-f64"))
+))]
+compile_error!("Feature `rapier` requires either `rapier-f32` or `rapier-f64` to be enabled.");
+
+#[cfg(all(feature = "dim2", feature = "rapier-f32"))]
 pub extern crate rapier2d as rapier;
-#[cfg(feature = "dim3")]
+#[cfg(all(feature = "dim2", feature = "rapier-f64"))]
+pub extern crate rapier2d_f64 as rapier;
+#[cfg(all(feature = "dim3", feature = "rapier-f32"))]
 pub extern crate rapier3d as rapier;
+#[cfg(all(feature = "dim3", feature = "rapier-f64"))]
+pub extern crate rapier3d_f64 as rapier;
 
 use slang_hal::re_exports::include_dir;
 
@@ -108,8 +124,43 @@ pub fn register_shaders(compiler: &mut SlangCompiler) {
 /// Re-exports Rapier's math types and defines dimension-specific type aliases
 /// for GPU simulation and angular inertia calculations.
 pub mod math {
-    /// Re-export all mathematical types from Rapier (vectors, matrices, etc.)
-    pub use rapier::math::*;
+    /// Scalar type used by the simulation.
+    pub type Real = f32;
+
+    /// Spatial point type.
+    #[cfg(feature = "dim2")]
+    pub type Point = glam::Vec2;
+    /// Spatial point type.
+    #[cfg(feature = "dim3")]
+    pub type Point = glam::Vec3;
+
+    /// Spatial vector type.
+    #[cfg(feature = "dim2")]
+    pub type Vector = glam::Vec2;
+    /// Spatial vector type.
+    #[cfg(feature = "dim3")]
+    pub type Vector = glam::Vec3;
+
+    /// Square matrix type.
+    #[cfg(feature = "dim2")]
+    pub type Matrix = glam::Mat2;
+    /// Square matrix type.
+    #[cfg(feature = "dim3")]
+    pub type Matrix = glam::Mat3;
+
+    /// Angular vector type.
+    #[cfg(feature = "dim2")]
+    pub type AngVector = f32;
+    /// Angular vector type.
+    #[cfg(feature = "dim3")]
+    pub type AngVector = glam::Vec3;
+
+    /// Spatial dimension.
+    #[cfg(feature = "dim2")]
+    pub const DIM: usize = 2;
+    /// Spatial dimension.
+    #[cfg(feature = "dim3")]
+    pub const DIM: usize = 3;
 
     /// GPU similarity transformation for 2D simulations (translation + rotation).
     #[cfg(feature = "dim2")]
@@ -124,6 +175,68 @@ pub mod math {
     /// Angular inertia type for 3D simulations (3x3 matrix).
     #[cfg(feature = "dim3")]
     pub type AngularInertia = glam::Mat3;
+
+    /// Conversions from Rapier's math types to the simulation's `f32` types.
+    ///
+    /// This is needed to support using Rapier with `f64` precision.
+    #[cfg(feature = "rapier")]
+    pub use rapier_convert::*;
+
+    #[cfg(feature = "rapier")]
+    mod rapier_convert {
+        use super::{Matrix, Vector};
+
+        /// Converts a Rapier scalar into the simulation's `f32` scalar.
+        #[cfg(feature = "rapier-f32")]
+        #[inline]
+        pub fn real(x: rapier::math::Real) -> f32 {
+            x
+        }
+        /// Converts a Rapier scalar into the simulation's `f32` scalar.
+        #[cfg(feature = "rapier-f64")]
+        #[inline]
+        pub fn real(x: rapier::math::Real) -> f32 {
+            x as f32
+        }
+
+        /// Converts a Rapier vector into the simulation's `f32` vector.
+        #[cfg(feature = "rapier-f32")]
+        #[inline]
+        pub fn vector(v: rapier::math::Vector) -> Vector {
+            v
+        }
+        /// Converts a Rapier vector into the simulation's `f32` vector.
+        #[cfg(all(feature = "rapier-f64", feature = "dim2"))]
+        #[inline]
+        pub fn vector(v: rapier::math::Vector) -> Vector {
+            v.as_vec2()
+        }
+        /// Converts a Rapier vector into the simulation's `f32` vector.
+        #[cfg(all(feature = "rapier-f64", feature = "dim3"))]
+        #[inline]
+        pub fn vector(v: rapier::math::Vector) -> Vector {
+            v.as_vec3()
+        }
+
+        /// Converts a Rapier matrix into the simulation's `f32` matrix.
+        #[cfg(feature = "rapier-f32")]
+        #[inline]
+        pub fn matrix(m: rapier::math::Matrix) -> Matrix {
+            m
+        }
+        /// Converts a Rapier matrix into the simulation's `f32` matrix.
+        #[cfg(all(feature = "rapier-f64", feature = "dim2"))]
+        #[inline]
+        pub fn matrix(m: rapier::math::Matrix) -> Matrix {
+            m.as_mat2()
+        }
+        /// Converts a Rapier matrix into the simulation's `f32` matrix.
+        #[cfg(all(feature = "rapier-f64", feature = "dim3"))]
+        #[inline]
+        pub fn matrix(m: rapier::math::Matrix) -> Matrix {
+            m.as_mat3()
+        }
+    }
 }
 
 /// Re-exports of commonly used dependencies for convenience.

src/pipeline.rs

@@ -8,6 +8,7 @@ use crate::grid::prefix_sum::{PrefixSumWorkspace, WgPrefixSum};
 use crate::grid::sort::WgSort;
 use crate::math::{GpuSim, Vector};
 use crate::rbd::dynamics::GpuBodySet;
+#[cfg(feature = "rapier")]
 use crate::rbd::dynamics::body::{BodyCoupling, BodyCouplingEntry};
 use crate::solver::{
     GpuBoundaryCondition, GpuImpulses, GpuMaterials, GpuParticleModelData, GpuParticles,
@@ -18,7 +19,9 @@ use crate::solver::{
 // The CDF kernel wrappers read the gated `Node.cdf`, so they only exist under the `cpic` feature.
 #[cfg(feature = "cpic")]
 use crate::solver::{WgG2PCdf, WgGridUpdateCdf, WgP2GCdf};
+#[cfg(feature = "rapier")]
 use rapier::dynamics::RigidBodySet;
+#[cfg(feature = "rapier")]
 use rapier::geometry::{ColliderHandle, ColliderSet};
 use slang_hal::backend::{Backend, Encoder, GpuTimestamps};
 use slang_hal::{BufferUsages, Shader, SlangCompiler};
@@ -275,6 +278,7 @@ pub struct MpmData<B: Backend, GpuModel: GpuParticleModelData> {
     /// Staging buffer for reading the timestep bound estimate.
     pub timestep_bounds_staging: GpuScalar<GpuTimestepBounds, B>,
     prefix_sum: PrefixSumWorkspace<B>,
+    #[cfg(feature = "rapier")]
     coupling: Vec<BodyCouplingEntry>,
 }
 
@@ -308,6 +312,7 @@ impl<B: Backend, GpuModel: GpuParticleModelData> MpmData<B, GpuModel> {
     /// # Returns
     ///
     /// GPU-resident simulation state ready for stepping.
+    #[cfg(feature = "rapier")]
     pub fn new(
         backend: &B,
         params: SimulationParams,
@@ -371,6 +376,7 @@ impl<B: Backend, GpuModel: GpuParticleModelData> MpmData<B, GpuModel> {
     /// # Returns
     ///
     /// GPU-resident simulation state ready for stepping.
+    #[cfg(feature = "rapier")]
     pub fn with_select_coupling(
         backend: &B,
         params: SimulationParams,
@@ -429,10 +435,78 @@ impl<B: Backend, GpuModel: GpuParticleModelData> MpmData<B, GpuModel> {
         })
     }
 
+    /// Creates MPM simulation data from a pre-built GPU rigid-body set.
+    ///
+    /// # Arguments
+    ///
+    /// * `backend` - GPU backend for buffer allocation
+    /// * `params` - Global simulation parameters (gravity, timestep)
+    /// * `particles` - Initial CPU-side particle data to upload
+    /// * `bodies` - GPU rigid bodies coupled with the simulation
+    /// * `materials` - Boundary condition per body (must have the same length as `bodies`)
+    /// * `cell_width` - Spatial width of each grid cell
+    /// * `grid_capacity` - Maximum number of active grid cells
+    ///
+    /// # Returns
+    ///
+    /// GPU-resident simulation state ready for stepping.
+    pub fn with_bodies(
+        backend: &B,
+        params: SimulationParams,
+        particles: &[Particle<GpuModel::Model>],
+        bodies: GpuBodySet<B>,
+        materials: &[GpuBoundaryCondition],
+        cell_width: f32,
+        grid_capacity: u32,
+    ) -> Result<Self, B::Error> {
+        let body_materials = GpuMaterials::new(backend, materials)?;
+        let sim_params = GpuSimulationParams::new(backend, params)?;
+        let particles = GpuParticles::from_particles(backend, particles)?;
+        let rigid_particles = GpuRigidParticles::new(backend)?;
+        let grid = GpuGrid::with_capacity(backend, grid_capacity, cell_width)?;
+        let prefix_sum = PrefixSumWorkspace::with_capacity(backend, grid_capacity)?;
+        let impulses = GpuImpulses::new(backend)?;
+        let poses_staging = GpuVector::vector_uninit(
+            backend,
+            bodies.len(),
+            BufferUsages::COPY_DST | BufferUsages::MAP_READ,
+        )?;
+        let bounds = GpuTimestepBounds::new();
+        let timestep_bounds = GpuTensor::scalar(
+            backend,
+            bounds,
+            BufferUsages::STORAGE | BufferUsages::COPY_SRC,
+        )?;
+        let timestep_bounds_staging = GpuTensor::scalar(
+            backend,
+            bounds,
+            BufferUsages::COPY_DST | BufferUsages::MAP_READ,
+        )?;
+
+        Ok(Self {
+            sim_params,
+            particles,
+            gravity: params.gravity,
+            rigid_particles,
+            bodies,
+            body_materials,
+            impulses,
+            grid,
+            prefix_sum,
+            poses_staging,
+            #[cfg(feature = "rapier")]
+            coupling: Vec::new(),
+            timestep_bounds,
+            timestep_bounds_staging,
+            base_dt: params.dt,
+        })
+    }
+
     /// Returns the list of rigid body coupling entries.
     ///
     /// Each entry specifies a collider-body pair that participates in MPM-rigid body
     /// interaction and the coupling mode.
+    #[cfg(feature = "rapier")]
     pub fn coupling(&self) -> &[BodyCouplingEntry] {
         &self.coupling
     }