diff --git a/CHANGELOG.md b/CHANGELOG.md
index 46f345d..d48d0a6 100644
--- a/CHANGELOG.md
+++ b/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 
diff --git a/crates/slosh2d/Cargo.toml b/crates/slosh2d/Cargo.toml
index f0d8aa3..36db2fc 100644
--- a/crates/slosh2d/Cargo.toml
+++ b/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"] }
diff --git a/crates/slosh3d/Cargo.toml b/crates/slosh3d/Cargo.toml
index 1cb8ad4..cbd6b57 100644
--- a/crates/slosh3d/Cargo.toml
+++ b/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"
diff --git a/crates/slosh_testbed2d/Cargo.toml b/crates/slosh_testbed2d/Cargo.toml
index dcd0a87..99a4fe2 100644
--- a/crates/slosh_testbed2d/Cargo.toml
+++ b/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"
 
diff --git a/crates/slosh_testbed3d/Cargo.toml b/crates/slosh_testbed3d/Cargo.toml
index 2d7f963..3c887c2 100644
--- a/crates/slosh_testbed3d/Cargo.toml
+++ b/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"
 
diff --git a/src/lib.rs b/src/lib.rs
index e978134..4a7eb09 100644
--- a/src/lib.rs
+++ b/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.
diff --git a/src/pipeline.rs b/src/pipeline.rs
index b71063b..2f704da 100644
--- a/src/pipeline.rs
+++ b/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
     }
diff --git a/src/rbd/dynamics/body.rs b/src/rbd/dynamics/body.rs
index df9be3e..3001120 100644
--- a/src/rbd/dynamics/body.rs
+++ b/src/rbd/dynamics/body.rs
@@ -2,8 +2,11 @@
 
 use crate::math::{AngVector, AngularInertia, GpuSim, Vector};
 use crate::rbd::shapes::{GpuShape, ShapeBuffers};
+#[cfg(feature = "rapier")]
 use rapier::geometry::ColliderHandle;
+#[cfg(feature = "rapier")]
 use rapier::prelude::MassProperties;
+#[cfg(feature = "rapier")]
 use rapier::{
     dynamics::{RigidBodyHandle, RigidBodySet},
     geometry::ColliderSet,
@@ -43,15 +46,17 @@ pub struct GpuMassProperties {
     pub com: Vector,
 }
 
+#[cfg(feature = "rapier")]
 impl From<MassProperties> for GpuMassProperties {
     fn from(props: MassProperties) -> Self {
+        use crate::math::{real, vector};
         GpuMassProperties {
             #[cfg(feature = "dim2")]
-            inv_inertia: props.inv_principal_inertia,
+            inv_inertia: real(props.inv_principal_inertia),
             #[cfg(feature = "dim3")]
-            inv_inertia: props.reconstruct_inverse_inertia_matrix(),
-            inv_mass: Vector::splat(props.inv_mass),
-            com: props.local_com,
+            inv_inertia: crate::math::matrix(props.reconstruct_inverse_inertia_matrix()),
+            inv_mass: Vector::splat(real(props.inv_mass)),
+            com: vector(props.local_com),
         }
     }
 }
@@ -139,6 +144,7 @@ pub enum BodyCoupling {
 ///
 /// Defines which Rapier rigid body and collider pair should be included in the
 /// GPU simulation and how they should be coupled.
+#[cfg(feature = "rapier")]
 #[derive(Copy, Clone, Debug, PartialEq, Eq)]
 pub struct BodyCouplingEntry {
     /// Handle to the Rapier rigid body
@@ -180,12 +186,15 @@ impl<B: Backend> GpuBodySet<B> {
     ///
     /// # Panics
     /// Panics if a collider has an unsupported shape type
+    #[cfg(feature = "rapier")]
     pub fn from_rapier(
         backend: &B,
         bodies: &RigidBodySet,
         colliders: &ColliderSet,
         coupling: &[BodyCouplingEntry],
     ) -> Result<Self, B::Error> {
+        use crate::math::vector;
+
         let mut shape_buffers = ShapeBuffers::default();
         let mut gpu_bodies = vec![];
         let mut pt_collider_ids = vec![];
@@ -205,14 +214,19 @@ impl<B: Backend> GpuBodySet<B> {
             let two_ways_coupling = rb.is_dynamic() && coupling.mode == BodyCoupling::TwoWays;
             let desc = BodyDesc {
                 vel: GpuVelocity {
-                    linear: rb.linvel(),
-                    #[allow(clippy::clone_on_copy)] // Needed for 2D/3D switch.
-                    angular: rb.angvel().clone(),
+                    linear: vector(rb.linvel()),
+                    #[cfg(feature = "dim2")]
+                    angular: crate::math::real(rb.angvel()),
+                    #[cfg(feature = "dim3")]
+                    angular: vector(rb.angvel()),
                 },
                 #[cfg(feature = "dim2")]
-                pose: GpuSim::from(rapier::na::Isometry2::from(*rb.position())),
+                pose: GpuSim::from(rapier::na::Isometry2::from(*rb.position()).cast::<f32>()),
                 #[cfg(feature = "dim3")]
-                pose: GpuSim::from_isometry((*rb.position()).into(), 1.0),
+                pose: GpuSim::from_isometry(
+                    rapier::na::Isometry3::from(*rb.position()).cast::<f32>(),
+                    1.0,
+                ),
                 shape,
                 local_mprops: if two_ways_coupling {
                     rb.mass_properties().local_mprops.into()
diff --git a/src/rbd/dynamics/mod.rs b/src/rbd/dynamics/mod.rs
index 35771a4..d35b56e 100644
--- a/src/rbd/dynamics/mod.rs
+++ b/src/rbd/dynamics/mod.rs
@@ -1,8 +1,8 @@
 //! Rigid-body dynamics (forces, velocities, etc.)
 
-pub use body::{
-    BodyCoupling, BodyCouplingEntry, BodyDesc, GpuBodySet, GpuForce, GpuMassProperties, GpuVelocity,
-};
+#[cfg(feature = "rapier")]
+pub use body::BodyCouplingEntry;
+pub use body::{BodyCoupling, BodyDesc, GpuBodySet, GpuForce, GpuMassProperties, GpuVelocity};
 
 /// Rigid body definitions and GPU body set management.
 pub mod body;
diff --git a/src/rbd/shapes.rs b/src/rbd/shapes.rs
index 126239d..3555437 100644
--- a/src/rbd/shapes.rs
+++ b/src/rbd/shapes.rs
@@ -5,6 +5,7 @@
 //! shapes to GPU-friendly formats with vertex buffers.
 
 use glam::Vec4;
+#[cfg(feature = "rapier")]
 use rapier::geometry::{Shape, ShapeType, TypedShape};
 
 use crate::math::Vector;
@@ -211,18 +212,23 @@ impl GpuShape {
     ///
     /// # Returns
     /// `Some(GpuShape)` if the shape type is supported, `None` otherwise
+    #[cfg(feature = "rapier")]
     pub fn from_parry(shape: &(impl Shape + ?Sized), buffers: &mut ShapeBuffers) -> Option<Self> {
+        use crate::math::{real, vector};
+
         match shape.as_typed_shape() {
-            TypedShape::Ball(shape) => Some(Self::ball(shape.radius)),
-            TypedShape::Cuboid(shape) => Some(Self::cuboid(shape.half_extents)),
+            TypedShape::Ball(shape) => Some(Self::ball(real(shape.radius))),
+            TypedShape::Cuboid(shape) => Some(Self::cuboid(vector(shape.half_extents))),
             TypedShape::Capsule(shape) => Some(Self::capsule(
-                shape.segment.a,
-                shape.segment.b,
-                shape.radius,
+                vector(shape.segment.a),
+                vector(shape.segment.b),
+                real(shape.radius),
             )),
             TypedShape::Polyline(shape) => {
                 let base_id = buffers.vertices.len();
-                buffers.vertices.extend_from_slice(shape.vertices());
+                buffers
+                    .vertices
+                    .extend(shape.vertices().iter().map(|v| vector(*v)));
                 Some(Self::polyline([
                     base_id as u32,
                     buffers.vertices.len() as u32,
@@ -230,7 +236,9 @@ impl GpuShape {
             }
             TypedShape::TriMesh(shape) => {
                 let base_id = buffers.vertices.len() as u32;
-                buffers.vertices.extend_from_slice(shape.vertices());
+                buffers
+                    .vertices
+                    .extend(shape.vertices().iter().map(|v| vector(*v)));
                 let gpu_trimesh = crate::trimesh::convert_trimesh_to_gpu(
                     shape,
                     &mut buffers.collision_vertices,
@@ -243,7 +251,7 @@ impl GpuShape {
             TypedShape::HeightField(shape) => {
                 let base_id = buffers.vertices.len();
                 let (vtx, _) = shape.to_polyline();
-                buffers.vertices.extend_from_slice(&vtx);
+                buffers.vertices.extend(vtx.iter().map(|v| vector(*v)));
                 Some(Self::polyline([
                     base_id as u32,
                     buffers.vertices.len() as u32,
@@ -253,7 +261,7 @@ impl GpuShape {
             TypedShape::HeightField(shape) => {
                 let (vtx, idx) = shape.to_trimesh();
                 let base_id = buffers.vertices.len() as u32;
-                buffers.vertices.extend_from_slice(&vtx);
+                buffers.vertices.extend(vtx.iter().map(|v| vector(*v)));
                 let trimesh = rapier::geometry::TriMesh::with_flags(
                     vtx,
                     idx,
@@ -269,9 +277,13 @@ impl GpuShape {
                 Some(Self::trimesh(&gpu_trimesh, base_id))
             }
             #[cfg(feature = "dim3")]
-            TypedShape::Cone(shape) => Some(Self::cone(shape.half_height, shape.radius)),
+            TypedShape::Cone(shape) => {
+                Some(Self::cone(real(shape.half_height), real(shape.radius)))
+            }
             #[cfg(feature = "dim3")]
-            TypedShape::Cylinder(shape) => Some(Self::cylinder(shape.half_height, shape.radius)),
+            TypedShape::Cylinder(shape) => {
+                Some(Self::cylinder(real(shape.half_height), real(shape.radius)))
+            }
             _ => None,
         }
     }
@@ -285,6 +297,7 @@ impl GpuShape {
     ///
     /// # Panics
     /// Panics if the stored type tag is invalid
+    #[cfg(feature = "rapier")]
     pub fn shape_type(&self) -> ShapeType {
         let tag = self.a.w.to_bits();
 
@@ -309,6 +322,7 @@ impl GpuShape {
     ///
     /// # Panics
     /// Panics if this shape is not a polyline
+    #[cfg(feature = "rapier")]
     pub fn polyline_rngs(&self) -> [u32; 2] {
         assert!(self.shape_type() == ShapeType::Polyline);
         [self.a.x.to_bits(), self.a.y.to_bits()]
@@ -318,6 +332,7 @@ impl GpuShape {
     ///
     /// # Panics
     /// Panics if this shape is not a triangle mesh
+    #[cfg(feature = "rapier")]
     pub fn trimesh_meta(&self) -> crate::trimesh::GpuTriMesh {
         assert!(self.shape_type() == ShapeType::TriMesh);
         crate::trimesh::GpuTriMesh {
@@ -333,6 +348,7 @@ impl GpuShape {
     ///
     /// # Panics
     /// Panics if this shape is not a triangle mesh
+    #[cfg(feature = "rapier")]
     pub fn trimesh_vertex_base_id(&self) -> u32 {
         assert!(self.shape_type() == ShapeType::TriMesh);
         self.b.z.to_bits()
diff --git a/src/sampling/mod.rs b/src/sampling/mod.rs
index 4032e2e..268a63f 100644
--- a/src/sampling/mod.rs
+++ b/src/sampling/mod.rs
@@ -4,12 +4,25 @@
 //! MPM-rigid body coupling. In 2D, samples polyline edges; in 3D, samples
 //! triangle mesh surfaces.
 
-#[cfg(feature = "dim2")]
-pub use sample_polyline::*;
-#[cfg(feature = "dim3")]
-pub use sample_trimesh::*;
+#[cfg(feature = "rapier")]
+mod rapier;
+#[cfg(feature = "rapier")]
+pub use rapier::*;
+
+use encase::ShaderType;
 
 #[cfg(feature = "dim2")]
-mod sample_polyline;
+#[derive(Copy, Clone, Debug, ShaderType)]
+#[repr(C)]
+pub struct GpuSampleIds {
+    pub segment: glam::UVec2,
+    pub collider: u32,
+}
+
 #[cfg(feature = "dim3")]
-mod sample_trimesh;
+#[derive(Copy, Clone, Debug, ShaderType)]
+#[repr(C)]
+pub struct GpuSampleIds {
+    pub triangle: glam::UVec3,
+    pub collider: u32,
+}
diff --git a/src/sampling/rapier/mod.rs b/src/sampling/rapier/mod.rs
new file mode 100644
index 0000000..c1e21b6
--- /dev/null
+++ b/src/sampling/rapier/mod.rs
@@ -0,0 +1,16 @@
+#[cfg(feature = "dim2")]
+mod sample_polyline;
+#[cfg(feature = "dim3")]
+mod sample_trimesh;
+
+#[cfg(feature = "dim2")]
+pub use sample_polyline::*;
+#[cfg(feature = "dim3")]
+pub use sample_trimesh::*;
+
+#[derive(Copy, Clone, Debug)]
+pub struct SamplingParams {
+    pub base_vid: u32,
+    pub collider_id: u32,
+    pub sampling_step: f32,
+}
diff --git a/src/sampling/sample_polyline.rs b/src/sampling/rapier/sample_polyline.rs
similarity index 52%
rename from src/sampling/sample_polyline.rs
rename to src/sampling/rapier/sample_polyline.rs
index 45248d2..c301150 100644
--- a/src/sampling/sample_polyline.rs
+++ b/src/sampling/rapier/sample_polyline.rs
@@ -1,22 +1,10 @@
-use crate::math::Vector;
-use encase::ShaderType;
-use glam::UVec2;
-use rapier::geometry::{Polyline, Segment};
-
-#[derive(Copy, Clone, Debug, ShaderType)]
-#[repr(C)]
-pub struct GpuSampleIds {
-    pub segment: UVec2,
-    pub collider: u32,
-}
+use crate::{
+    math::{Vector, vector},
+    sampling::{GpuSampleIds, SamplingParams},
+};
 
-#[derive(Copy, Clone, Debug)]
-#[repr(C)]
-pub struct SamplingParams {
-    pub base_vid: u32,
-    pub collider_id: u32,
-    pub sampling_step: f32,
-}
+use glam::UVec2;
+use rapier::geometry::Polyline;
 
 #[derive(Default, Clone)]
 pub struct SamplingBuffers {
@@ -30,29 +18,30 @@ pub fn sample_polyline(
     buffers: &mut SamplingBuffers,
 ) {
     for seg_idx in polyline.indices() {
-        let seg = Segment::new(
-            polyline.vertices()[seg_idx[0] as usize],
-            polyline.vertices()[seg_idx[1] as usize],
-        );
+        let a = vector(polyline.vertices()[seg_idx[0] as usize]);
+        let b = vector(polyline.vertices()[seg_idx[1] as usize]);
         let sample_id = GpuSampleIds {
             segment: UVec2::new(params.base_vid + seg_idx[0], params.base_vid + seg_idx[1]),
             collider: params.collider_id,
         };
-        buffers.samples.push(seg.a);
+        buffers.samples.push(a);
         buffers.samples_ids.push(sample_id);
 
-        if let Some(dir) = seg.direction() {
+        let ab = b - a;
+        let length = ab.length();
+        if length > 0.0 {
+            let dir = ab / length;
             for i in 0.. {
                 let shift = (i as f32) * params.sampling_step;
-                if shift > seg.length() {
+                if shift > length {
                     break;
                 }
 
-                buffers.samples.push(seg.a + dir * shift);
+                buffers.samples.push(a + dir * shift);
                 buffers.samples_ids.push(sample_id);
             }
 
-            buffers.samples.push(seg.b);
+            buffers.samples.push(b);
             buffers.samples_ids.push(sample_id);
         }
     }
diff --git a/src/sampling/sample_trimesh.rs b/src/sampling/rapier/sample_trimesh.rs
similarity index 87%
rename from src/sampling/sample_trimesh.rs
rename to src/sampling/rapier/sample_trimesh.rs
index 165d493..21c3b89 100644
--- a/src/sampling/sample_trimesh.rs
+++ b/src/sampling/rapier/sample_trimesh.rs
@@ -1,41 +1,44 @@
-use crate::math::Vector;
-use encase::ShaderType;
+use crate::{
+    math::{Vector, vector},
+    sampling::{GpuSampleIds, SamplingParams},
+};
+
 use glam::UVec3;
-use rapier::geometry::{Segment, TriMesh, Triangle};
+use rapier::geometry::TriMesh;
 use std::collections::HashSet;
 
+#[derive(Default, Clone)]
+pub struct SamplingBuffers {
+    pub local_samples: Vec<Vector>,
+    pub samples: Vec<Vector>,
+    pub samples_ids: Vec<GpuSampleIds>,
+}
+
 // Epsilon used as a length threshold in various steps of the sampling. In particular, this avoids
 // degenerate geometries from generating invalid samples.
 const EPS: f32 = 1.0e-5;
 
-pub struct TriangleSample {
-    pub triangle_id: u32,
-    pub point: Vector,
-}
-
-#[derive(Copy, Clone, Debug, ShaderType)]
-#[repr(C)]
-pub struct GpuSampleIds {
-    pub triangle: UVec3,
-    pub collider: u32,
+#[derive(Copy, Clone)]
+struct Segment {
+    a: Vector,
+    b: Vector,
 }
 
-#[derive(Copy, Clone, Debug)]
-pub struct SamplingParams {
-    pub base_vid: u32,
-    pub collider_id: u32,
-    pub sampling_step: f32,
+#[derive(Copy, Clone)]
+struct Triangle {
+    a: Vector,
+    b: Vector,
+    c: Vector,
 }
 
-#[derive(Default, Clone)]
-pub struct SamplingBuffers {
-    pub local_samples: Vec<Vector>,
-    pub samples: Vec<Vector>,
-    pub samples_ids: Vec<GpuSampleIds>,
+pub struct TriangleSample {
+    pub triangle_id: u32,
+    pub point: Vector,
 }
 
 pub fn sample_trimesh(trimesh: &TriMesh, params: &SamplingParams, buffers: &mut SamplingBuffers) {
-    let samples = sample_mesh(trimesh.vertices(), trimesh.indices(), params.sampling_step);
+    let vertices: Vec<Vector> = trimesh.vertices().iter().map(|v| vector(*v)).collect();
+    let samples = sample_mesh(&vertices, trimesh.indices(), params.sampling_step);
 
     for sample in samples {
         let tri_idx = trimesh.indices()[sample.triangle_id as usize];
@@ -80,25 +83,34 @@ pub fn sample_mesh(
     };
 
     for (tri_id, idx) in indices.iter().enumerate() {
-        let tri = Triangle::new(
-            vertices[idx[0] as usize],
-            vertices[idx[1] as usize],
-            vertices[idx[2] as usize],
-        );
+        let tri = Triangle {
+            a: vertices[idx[0] as usize],
+            b: vertices[idx[1] as usize],
+            c: vertices[idx[2] as usize],
+        };
         sample_triangle(tri, &mut samples, xy_spacing, tri_id as u32);
 
         if seg_needs_sampling(idx[0], idx[1]) {
-            let seg = Segment::new(vertices[idx[0] as usize], vertices[idx[1] as usize]);
+            let seg = Segment {
+                a: vertices[idx[0] as usize],
+                b: vertices[idx[1] as usize],
+            };
             sample_edge(seg, &mut samples, xy_spacing, tri_id as u32);
         }
 
         if seg_needs_sampling(idx[1], idx[2]) {
-            let seg = Segment::new(vertices[idx[1] as usize], vertices[idx[2] as usize]);
+            let seg = Segment {
+                a: vertices[idx[1] as usize],
+                b: vertices[idx[2] as usize],
+            };
             sample_edge(seg, &mut samples, xy_spacing, tri_id as u32);
         }
 
         if seg_needs_sampling(idx[2], idx[0]) {
-            let seg = Segment::new(vertices[idx[2] as usize], vertices[idx[0] as usize]);
+            let seg = Segment {
+                a: vertices[idx[2] as usize],
+                b: vertices[idx[0] as usize],
+            };
             sample_edge(seg, &mut samples, xy_spacing, tri_id as u32);
         }
     }
@@ -111,7 +123,7 @@ pub fn sample_mesh(
 ///
 /// The returned samples will not contain `edge.a`. It might contain `edge.b` (but it is unlikely)
 /// if it aligns exactly with the internal sampling spacing.
-pub fn sample_edge(
+fn sample_edge(
     edge: Segment,
     samples: &mut Vec<TriangleSample>,
     xy_spacing: f32,
@@ -147,7 +159,7 @@ pub fn sample_edge(
 ///
 /// Because this does not attempt to sample the edges of the triangles, small or thin triangles
 /// might not result in any samples. Edges should be sampled separately with [`sample_edge`].
-pub fn sample_triangle(
+fn sample_triangle(
     triangle: Triangle,
     samples: &mut Vec<TriangleSample>,
     xy_spacing: f32,
diff --git a/src/solver/particle.rs b/src/solver/particle.rs
index 952be49..35c278b 100644
--- a/src/solver/particle.rs
+++ b/src/solver/particle.rs
@@ -1,12 +1,17 @@
 use crate::math::{Matrix, Vector};
+#[cfg(feature = "rapier")]
 use crate::rbd::dynamics::GpuBodySet;
+#[cfg(feature = "rapier")]
 use crate::rbd::dynamics::body::BodyCouplingEntry;
-use crate::rbd::shapes::ShapeBuffers;
+#[cfg(feature = "rapier")]
 use crate::sampling;
-use crate::sampling::{GpuSampleIds, SamplingBuffers, SamplingParams};
+use crate::sampling::GpuSampleIds;
+#[cfg(feature = "rapier")]
+use crate::sampling::{SamplingBuffers, SamplingParams};
 use crate::solver::particle_model::GpuParticleModelData;
 use bytemuck::{Pod, Zeroable};
 use encase::ShaderType;
+#[cfg(feature = "rapier")]
 use rapier::geometry::ColliderSet;
 use slang_hal::{BufferUsages, backend::Backend};
 use std::ops::RangeBounds;
@@ -271,16 +276,25 @@ pub struct GpuRigidParticles<B: Backend> {
 impl<B: Backend> GpuRigidParticles<B> {
     /// Creates an empty set of rigid particles.
     pub fn new(backend: &B) -> Result<Self, B::Error> {
-        Self::from_rapier(
-            backend,
-            &ColliderSet::default(),
-            &GpuBodySet::new(backend, &[], &[], &ShapeBuffers::default())?,
-            &[],
-            1.0,
-        )
+        Ok(Self {
+            local_sample_points: GpuTensor::vector_uninit_encased(
+                backend,
+                1,
+                BufferUsages::STORAGE,
+            )?,
+            sample_points: GpuTensor::vector_uninit_encased(backend, 1, BufferUsages::STORAGE)?,
+            node_linked_lists: GpuTensor::vector_uninit(backend, 1, BufferUsages::STORAGE)?,
+            sample_ids: GpuTensor::vector_uninit_encased(backend, 1, BufferUsages::STORAGE)?,
+            rigid_particle_needs_block: GpuTensor::vector_uninit(
+                backend,
+                1,
+                BufferUsages::STORAGE,
+            )?,
+        })
     }
 
     /// Samples particles from Rapier collider surfaces for MPM coupling.
+    #[cfg(feature = "rapier")]
     pub fn from_rapier(
         backend: &B,
         colliders: &ColliderSet,
diff --git a/src/trimesh.rs b/src/trimesh.rs
index a6ffc2b..fb3af98 100644
--- a/src/trimesh.rs
+++ b/src/trimesh.rs
@@ -1,9 +1,8 @@
 // TODO: move this to rbd?
 
-use crate::math::Vector;
 use encase::ShaderType;
-use rapier::geometry::TriMesh;
-use rapier::prelude::DIM;
+#[cfg(feature = "rapier")]
+use {crate::math::Vector, rapier::geometry::TriMesh, rapier::prelude::DIM};
 
 #[derive(Copy, Clone, ShaderType)]
 pub struct GpuTriMesh {
@@ -19,11 +18,14 @@ pub struct GpuTriMesh {
     pub num_vertices: u32,
 }
 
+#[cfg(feature = "rapier")]
 pub fn convert_trimesh_to_gpu(
     shape: &TriMesh,
     vertices: &mut Vec<Vector>,
     indices: &mut Vec<u32>,
 ) -> GpuTriMesh {
+    use crate::math::vector;
+
     let bvh_vtx_root_id = vertices.len();
     let bvh_idx_root_id = indices.len();
     // Append the BVH data to the vertex/index buffers.
@@ -57,7 +59,10 @@ pub fn convert_trimesh_to_gpu(
         .map(|tri| {
             let aabb = tri.local_aabb();
             BvhObject {
-                aabb: bvh::aabb::Aabb::with_bounds(aabb.mins.into(), aabb.maxs.into()),
+                aabb: bvh::aabb::Aabb::with_bounds(
+                    vector(aabb.mins).into(),
+                    vector(aabb.maxs).into(),
+                ),
                 node_index: 0,
             }
         })
@@ -83,14 +88,10 @@ pub fn convert_trimesh_to_gpu(
         let pn = shape
             .pseudo_normals()
             .expect("trimeshes without pseudo-normals are not supported");
-        vertices.extend_from_slice(shape.vertices());
-        vertices.extend_from_slice(&pn.vertices_pseudo_normal);
+        vertices.extend(shape.vertices().iter().map(|v| vector(*v)));
+        vertices.extend(pn.vertices_pseudo_normal.iter().map(|v| vector(*v)));
         assert_eq!(shape.vertices().len(), pn.vertices_pseudo_normal.len());
-        vertices.extend(
-            pn.edges_pseudo_normal
-                .iter()
-                .flat_map(|n| n.map(Vector::from)),
-        );
+        vertices.extend(pn.edges_pseudo_normal.iter().flat_map(|n| n.map(vector)));
     }
     indices.extend(shape.indices().iter().flat_map(|tri| tri.iter().copied()));
     GpuTriMesh {