diff --git a/CHANGELOG.md b/CHANGELOG.md
index 8ea2906..46f345d 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,4 +1,5 @@
 # Unreleased
+- 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 
   kernels now clamp the dispatch and grid-stride over the particles.
diff --git a/Cargo.toml b/Cargo.toml
index 6b8c2c0..9276d02 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -4,11 +4,16 @@ members = [
     "crates/slosh_testbed3d",
     "crates/slosh2d",
     "crates/slosh3d",
-#    "validation_tests",
+    # "validation_tests",
 ]
 resolver = "2"
 
 [workspace.dependencies]
+glam = { version = "0.30", default-features = false, features = [
+    "std",
+    "encase",
+    "bytemuck",
+] }
 nalgebra = { version = "0.34", features = ["convert-bytemuck"] }
 wgpu = { version = "27", features = ["naga-ir"] }
 encase = "0.12"
diff --git a/crates/slosh2d/Cargo.toml b/crates/slosh2d/Cargo.toml
index 2c0523d..b73f94e 100644
--- a/crates/slosh2d/Cargo.toml
+++ b/crates/slosh2d/Cargo.toml
@@ -19,12 +19,12 @@ required-features = ["dim2"]
 default = ["dim2"]
 dim2 = []
 
-comptime = [ "slosh_testbed2d/comptime", "stensor/comptime" ]
-runtime = [ "slosh_testbed2d/runtime", "stensor/runtime" ]
-webgpu = [ "slosh_testbed2d/webgpu", "stensor/webgpu" ]
+comptime = ["slosh_testbed2d/comptime", "stensor/comptime"]
+runtime = ["slosh_testbed2d/runtime", "stensor/runtime"]
+webgpu = ["slosh_testbed2d/webgpu", "stensor/webgpu"]
 
 [dependencies]
-nalgebra = { workspace = true }
+glam = { workspace = true }
 wgpu = { workspace = true }
 encase = { workspace = true }
 slang-hal = { workspace = true }
@@ -37,7 +37,7 @@ static_assertions = { workspace = true }
 bvh = "0.12.0"
 
 # TODO: make rapier optional?
-rapier2d = "0.31"
+rapier2d = "0.32"
 
 # For wasm?
 getrandom = { version = "0.3.1", features = ["wasm_js"] }
@@ -45,12 +45,13 @@ uuid = { version = "1", features = ["js"] }
 
 
 [dev-dependencies]
-nalgebra = { version = "0.34", features = ["rand"] }
 futures-test = "0.3"
 serial_test = "3"
 approx = "0.5"
 async-std = { version = "1", features = ["attributes"] }
-slosh_testbed2d = { path = "../slosh_testbed2d", features = ["dim2"], default-features = false }
+slosh_testbed2d = { path = "../slosh_testbed2d", features = [
+    "dim2",
+], default-features = false }
 kiss3d = { version = "0.37.0" }
 
 [build-dependencies]
diff --git a/crates/slosh2d/examples/centilever_beam2.rs b/crates/slosh2d/examples/centilever_beam2.rs
index cebf182..11c99a2 100644
--- a/crates/slosh2d/examples/centilever_beam2.rs
+++ b/crates/slosh2d/examples/centilever_beam2.rs
@@ -1,6 +1,6 @@
 use slosh_testbed2d::{RapierData, slosh};
 
-use nalgebra::{point, vector};
+use glam::vec2;
 use rapier2d::prelude::{ColliderBuilder, RigidBodyBuilder};
 use slang_hal::backend::WebGpu;
 use slosh::{
@@ -33,7 +33,7 @@ pub fn beam_demo(backend: &WebGpu, app_state: &mut AppState) -> PhysicsContext {
     let mut particles = vec![];
     for i in 0..ni {
         for j in 0..nj {
-            let position = point![i as f32, j as f32] * diameter;
+            let position = vec2(i as f32, j as f32) * diameter;
             let density = 1000.0;
             let radius = diameter / 2.0;
             let model = ParticleModel::elastic_neo_hookean(young_modulus, poisson_ratio);
@@ -48,13 +48,13 @@ pub fn beam_demo(backend: &WebGpu, app_state: &mut AppState) -> PhysicsContext {
     };
 
     let params = SimulationParams {
-        gravity: vector![0.0, -9.81] * app_state.gravity_factor,
+        gravity: vec2(0.0, -9.81) * app_state.gravity_factor,
         dt: 1.0 / 60.0,
         padding: 0.0,
     };
 
     let rb = RigidBodyBuilder::fixed()
-        .translation(vector![0.0, height / 2.0])
+        .translation(vec2(0.0, height / 2.0))
         .build();
     let rb_handle = rapier_data.bodies.insert(rb);
     let co = ColliderBuilder::cuboid(fixed_part, height);
diff --git a/crates/slosh2d/examples/elastic_cut2.rs b/crates/slosh2d/examples/elastic_cut2.rs
index ae23d77..6a0f238 100644
--- a/crates/slosh2d/examples/elastic_cut2.rs
+++ b/crates/slosh2d/examples/elastic_cut2.rs
@@ -1,6 +1,6 @@
 use slosh_testbed2d::{RapierData, slosh};
 
-use nalgebra::{Vector2, point, vector};
+use glam::vec2;
 use rapier2d::prelude::{ColliderBuilder, RigidBodyBuilder};
 use slang_hal::backend::WebGpu;
 use slosh::{
@@ -25,7 +25,7 @@ pub fn elastic_cut_demo(backend: &WebGpu, app_state: &mut AppState) -> PhysicsCo
     for i in 0..700 {
         for j in 0..700 {
             let position =
-                point![i as f32 + 0.5, j as f32 + 0.5] * cell_width / 2.0 + Vector2::y() * offset_y;
+                vec2(i as f32 + 0.5, j as f32 + 0.5) * cell_width / 2.0 + vec2(0.0, offset_y);
 
             let density = 1000.0;
             let radius = cell_width / 4.0;
@@ -40,7 +40,7 @@ pub fn elastic_cut_demo(backend: &WebGpu, app_state: &mut AppState) -> PhysicsCo
     };
 
     let params = SimulationParams {
-        gravity: vector![0.0, -9.81] * app_state.gravity_factor,
+        gravity: vec2(0.0, -9.81) * app_state.gravity_factor,
         dt: 1.0 / 60.0,
         padding: 0.0,
     };
@@ -50,7 +50,7 @@ pub fn elastic_cut_demo(backend: &WebGpu, app_state: &mut AppState) -> PhysicsCo
     /*
      * Static platforms.
      */
-    let rb = RigidBodyBuilder::fixed().translation(vector![35.0, 20.0]);
+    let rb = RigidBodyBuilder::fixed().translation(vec2(35.0, 20.0));
     let rb_handle = rapier_data.bodies.insert(rb);
     let co = ColliderBuilder::cuboid(70.0, 1.0);
     rapier_data
@@ -60,12 +60,12 @@ pub fn elastic_cut_demo(backend: &WebGpu, app_state: &mut AppState) -> PhysicsCo
     let mut polyline = vec![];
     let subdivs = 100;
     let length = 84.0;
-    let start = point![35.0, 70.0] - vector![length / 2.0, 0.0];
+    let start = vec2(35.0 - length / 2.0, 70.0);
 
     for i in 0..=subdivs {
         let step = length / (subdivs as f32);
         let dx = i as f32 * step;
-        polyline.push(start + vector![dx, dx.sin()])
+        polyline.push(start + vec2(dx, dx.sin()))
     }
 
     let rb = RigidBodyBuilder::fixed();
@@ -80,8 +80,8 @@ pub fn elastic_cut_demo(backend: &WebGpu, app_state: &mut AppState) -> PhysicsCo
         let rb_handle = rapier_data.bodies.insert(rb);
         let co = ColliderBuilder::polyline(
             vec![
-                point![0.0 + k as f32 * 15.0, 20.0],
-                point![-10.0 + k as f32 * 15.0, 45.0],
+                vec2(0.0 + k as f32 * 15.0, 20.0),
+                vec2(-10.0 + k as f32 * 15.0, 45.0),
             ],
             None,
         )
diff --git a/crates/slosh2d/examples/elasticity2.rs b/crates/slosh2d/examples/elasticity2.rs
index e19174a..ee0d33a 100644
--- a/crates/slosh2d/examples/elasticity2.rs
+++ b/crates/slosh2d/examples/elasticity2.rs
@@ -1,6 +1,6 @@
 use slosh_testbed2d::{RapierData, slosh};
 
-use nalgebra::{Vector2, point, vector};
+use glam::vec2;
 use rapier2d::prelude::{ColliderBuilder, RigidBodyBuilder};
 use slang_hal::backend::WebGpu;
 use slosh::{
@@ -25,8 +25,8 @@ pub fn elasticity_demo(backend: &WebGpu, app_state: &mut AppState) -> PhysicsCon
     for i in 0..700 {
         for j in 0..700 {
             let position =
-                point![i as f32 + 0.5 + (i / 50) as f32 * 2.0, j as f32 + 0.5] * cell_width / 2.0
-                    + Vector2::y() * offset_y;
+                vec2(i as f32 + 0.5 + (i / 50) as f32 * 2.0, j as f32 + 0.5) * cell_width / 2.0
+                    + vec2(0.0, offset_y);
             let density = 1000.0;
             let radius = cell_width / 4.0;
             let model = ParticleModel::elastic(5.0e6, 0.2);
@@ -40,12 +40,12 @@ pub fn elasticity_demo(backend: &WebGpu, app_state: &mut AppState) -> PhysicsCon
     };
 
     let params = SimulationParams {
-        gravity: vector![0.0, -9.81] * app_state.gravity_factor,
+        gravity: vec2(0.0, -9.81) * app_state.gravity_factor,
         dt: 1.0 / 60.0,
         padding: 0.0,
     };
 
-    let rb = RigidBodyBuilder::fixed().translation(vector![0.0, -1.0]);
+    let rb = RigidBodyBuilder::fixed().translation(vec2(0.0, -1.0));
     let rb_handle = rapier_data.bodies.insert(rb);
     let co = ColliderBuilder::cuboid(1000.0, 1.0);
     rapier_data
@@ -53,7 +53,7 @@ pub fn elasticity_demo(backend: &WebGpu, app_state: &mut AppState) -> PhysicsCon
         .insert_with_parent(co, rb_handle, &mut rapier_data.bodies);
 
     let rb = RigidBodyBuilder::fixed()
-        .translation(vector![-20.0, 0.0])
+        .translation(vec2(-20.0, 0.0))
         .rotation(0.5);
     let rb_handle = rapier_data.bodies.insert(rb);
     let co = ColliderBuilder::cuboid(1.0, 60.0);
@@ -62,7 +62,7 @@ pub fn elasticity_demo(backend: &WebGpu, app_state: &mut AppState) -> PhysicsCon
         .insert_with_parent(co, rb_handle, &mut rapier_data.bodies);
 
     let rb = RigidBodyBuilder::fixed()
-        .translation(vector![90.0, 0.0])
+        .translation(vec2(90.0, 0.0))
         .rotation(-0.5);
     let rb_handle = rapier_data.bodies.insert(rb);
     let co = ColliderBuilder::cuboid(1.0, 60.0);
diff --git a/crates/slosh2d/examples/sand2.rs b/crates/slosh2d/examples/sand2.rs
index c60cf2c..a209681 100644
--- a/crates/slosh2d/examples/sand2.rs
+++ b/crates/slosh2d/examples/sand2.rs
@@ -1,6 +1,6 @@
 use slosh_testbed2d::{RapierData, slosh};
 
-use nalgebra::{Vector2, point, vector};
+use glam::vec2;
 use rapier2d::prelude::{ColliderBuilder, RigidBodyBuilder};
 use slang_hal::backend::WebGpu;
 use slosh::{
@@ -25,7 +25,7 @@ pub fn sand_demo(backend: &WebGpu, app_state: &mut AppState) -> PhysicsContext {
     for i in 0..700 {
         for j in 0..700 {
             let position =
-                point![i as f32 + 0.5, j as f32 + 0.5] * cell_width / 2.0 + Vector2::y() * offset_y;
+                vec2(i as f32 + 0.5, j as f32 + 0.5) * cell_width / 2.0 + vec2(0.0, offset_y);
             let density = 1000.0;
             let radius = cell_width / 4.0;
             let young_modulus = 1.0e7;
@@ -43,7 +43,7 @@ pub fn sand_demo(backend: &WebGpu, app_state: &mut AppState) -> PhysicsContext {
     };
 
     let params = SimulationParams {
-        gravity: vector![0.0, -9.81] * app_state.gravity_factor,
+        gravity: vec2(0.0, -9.81) * app_state.gravity_factor,
         dt: 1.0 / 60.0,
         padding: 0.0,
     };
@@ -53,7 +53,7 @@ pub fn sand_demo(backend: &WebGpu, app_state: &mut AppState) -> PhysicsContext {
     /*
      * Static platforms.
      */
-    let rb = RigidBodyBuilder::fixed().translation(vector![35.0, -1.0]);
+    let rb = RigidBodyBuilder::fixed().translation(vec2(35.0, -1.0));
     let rb_handle = rapier_data.bodies.insert(rb);
     let co = ColliderBuilder::cuboid(42.0, 1.0);
     rapier_data
@@ -61,7 +61,7 @@ pub fn sand_demo(backend: &WebGpu, app_state: &mut AppState) -> PhysicsContext {
         .insert_with_parent(co, rb_handle, &mut rapier_data.bodies);
 
     let rb = RigidBodyBuilder::fixed()
-        .translation(vector![-25.0, 45.0])
+        .translation(vec2(-25.0, 45.0))
         .rotation(0.5);
     let rb_handle = rapier_data.bodies.insert(rb);
     let co = ColliderBuilder::cuboid(1.0, 52.0);
@@ -70,7 +70,7 @@ pub fn sand_demo(backend: &WebGpu, app_state: &mut AppState) -> PhysicsContext {
         .insert_with_parent(co, rb_handle, &mut rapier_data.bodies);
 
     let rb = RigidBodyBuilder::fixed()
-        .translation(vector![95.0, 45.0])
+        .translation(vec2(95.0, 45.0))
         .rotation(-0.5);
     let rb_handle = rapier_data.bodies.insert(rb);
     let co = ColliderBuilder::cuboid(1.0, 52.0);
@@ -82,7 +82,7 @@ pub fn sand_demo(backend: &WebGpu, app_state: &mut AppState) -> PhysicsContext {
      * Rotating platforms.
      */
     let rb = RigidBodyBuilder::kinematic_velocity_based()
-        .translation(vector![5.0, 35.0])
+        .translation(vec2(5.0, 35.0))
         .angvel(ANGVEL);
     let rb_handle = rapier_data.bodies.insert(rb);
     let co = ColliderBuilder::cuboid(1.0, 10.0);
@@ -91,7 +91,7 @@ pub fn sand_demo(backend: &WebGpu, app_state: &mut AppState) -> PhysicsContext {
         .insert_with_parent(co, rb_handle, &mut rapier_data.bodies);
 
     let rb = RigidBodyBuilder::kinematic_velocity_based()
-        .translation(vector![35.0, 35.0])
+        .translation(vec2(35.0, 35.0))
         .angvel(-ANGVEL);
     let rb_handle = rapier_data.bodies.insert(rb);
     let co = ColliderBuilder::cuboid(10.0, 1.0);
@@ -100,7 +100,7 @@ pub fn sand_demo(backend: &WebGpu, app_state: &mut AppState) -> PhysicsContext {
         .insert_with_parent(co, rb_handle, &mut rapier_data.bodies);
 
     let rb = RigidBodyBuilder::kinematic_velocity_based()
-        .translation(vector![65.0, 35.0])
+        .translation(vec2(65.0, 35.0))
         .angvel(ANGVEL);
     let rb_handle = rapier_data.bodies.insert(rb);
     let co = ColliderBuilder::cuboid(1.0, 10.0);
@@ -109,7 +109,7 @@ pub fn sand_demo(backend: &WebGpu, app_state: &mut AppState) -> PhysicsContext {
         .insert_with_parent(co, rb_handle, &mut rapier_data.bodies);
 
     let rb = RigidBodyBuilder::kinematic_velocity_based()
-        .translation(vector![20.0, 20.0])
+        .translation(vec2(20.0, 20.0))
         .angvel(-ANGVEL);
     let rb_handle = rapier_data.bodies.insert(rb);
     let co = ColliderBuilder::ball(5.0);
@@ -118,7 +118,7 @@ pub fn sand_demo(backend: &WebGpu, app_state: &mut AppState) -> PhysicsContext {
         .insert_with_parent(co, rb_handle, &mut rapier_data.bodies);
 
     let rb = RigidBodyBuilder::kinematic_velocity_based()
-        .translation(vector![50.0, 20.0])
+        .translation(vec2(50.0, 20.0))
         .angvel(-ANGVEL);
     let rb_handle = rapier_data.bodies.insert(rb);
     let co = ColliderBuilder::capsule_y(5.0, 3.0);
diff --git a/crates/slosh3d/Cargo.toml b/crates/slosh3d/Cargo.toml
index 7c0aeff..b0c0b2d 100644
--- a/crates/slosh3d/Cargo.toml
+++ b/crates/slosh3d/Cargo.toml
@@ -19,12 +19,12 @@ required-features = ["dim3"]
 default = ["dim3"]
 dim3 = []
 
-comptime = [ "stensor/comptime", "slosh_testbed3d/comptime" ]
-runtime = [ "stensor/runtime", "slosh_testbed3d/runtime" ]
-webgpu = [ "stensor/webgpu", "slosh_testbed3d/webgpu" ]
+comptime = ["stensor/comptime", "slosh_testbed3d/comptime"]
+runtime = ["stensor/runtime", "slosh_testbed3d/runtime"]
+webgpu = ["stensor/webgpu", "slosh_testbed3d/webgpu"]
 
 [dependencies]
-nalgebra = { workspace = true }
+glam = { workspace = true }
 wgpu = { workspace = true }
 encase = { workspace = true }
 bytemuck = { workspace = true }
@@ -37,10 +37,9 @@ bvh = "0.12.0"
 serde = "1"
 
 # TODO: make rapier optional?
-rapier3d = "0.31"
+rapier3d = "0.32"
 
 [dev-dependencies]
-nalgebra = { version = "0.34", features = ["rand"] }
 futures-test = "0.3"
 serial_test = "3"
 approx = "0.5"
diff --git a/crates/slosh3d/examples/centilever_beam3.rs b/crates/slosh3d/examples/centilever_beam3.rs
index c988e4c..bcba268 100644
--- a/crates/slosh3d/examples/centilever_beam3.rs
+++ b/crates/slosh3d/examples/centilever_beam3.rs
@@ -1,6 +1,6 @@
 use slosh_testbed3d::{PhysicsState, RapierData, slosh};
 
-use nalgebra::{point, vector};
+use glam::vec3;
 use rapier3d::prelude::{ColliderBuilder, RigidBodyBuilder};
 use slang_hal::backend::WebGpu;
 use slosh::{
@@ -34,7 +34,7 @@ pub fn beam_demo(backend: &WebGpu, app_state: &mut AppState) -> PhysicsContext {
     for i in 0..ni {
         for j in 0..njk {
             for k in 0..njk {
-                let position = point![i as f32, j as f32, k as f32] * diameter;
+                let position = vec3(i as f32, j as f32, k as f32) * diameter;
                 let density = 1000.0;
                 let radius = diameter / 2.0;
                 let model = ParticleModel::elastic_neo_hookean(young_modulus, poisson_ratio);
@@ -52,12 +52,12 @@ pub fn beam_demo(backend: &WebGpu, app_state: &mut AppState) -> PhysicsContext {
     };
 
     let params = SimulationParams {
-        gravity: vector![0.0, -9.81, 0.0] * app_state.gravity_factor,
+        gravity: vec3(0.0, -9.81, 0.0) * app_state.gravity_factor,
         dt: 1.0 / 60.0,
     };
 
     let rb = RigidBodyBuilder::fixed()
-        .translation(vector![0.0, height / 2.0, height / 2.0])
+        .translation(vec3(0.0, height / 2.0, height / 2.0))
         .build();
     let rb_handle = rapier_data.bodies.insert(rb);
     let co = ColliderBuilder::cuboid(fixed_part, height, height);
diff --git a/crates/slosh3d/examples/elastic_cut3.rs b/crates/slosh3d/examples/elastic_cut3.rs
index cd6e09e..71fe7f0 100644
--- a/crates/slosh3d/examples/elastic_cut3.rs
+++ b/crates/slosh3d/examples/elastic_cut3.rs
@@ -1,7 +1,8 @@
 use slosh_testbed3d::{RapierData, slosh};
 
-use nalgebra::{DMatrix, Isometry3, point, vector};
+use glam::{Quat, vec3};
 use rapier3d::geometry::HeightField;
+use rapier3d::parry::utils::Array2;
 use rapier3d::prelude::{ColliderBuilder, RigidBodyBuilder};
 use slang_hal::backend::WebGpu;
 use slosh::{
@@ -24,11 +25,11 @@ pub fn elastic_cut_demo(backend: &WebGpu, app_state: &mut AppState) -> PhysicsCo
     for i in 0..nxz {
         for j in 0..30 {
             for k in 0..nxz {
-                let position = point![
+                let position = vec3(
                     i as f32 + 0.5 - nxz as f32 / 2.0,
                     j as f32 + 0.5 + 60.0,
-                    k as f32 + 0.5 - nxz as f32 / 2.0
-                ] * cell_width
+                    k as f32 + 0.5 - nxz as f32 / 2.0,
+                ) * cell_width
                     / 2.0;
                 let density = 2700.0;
                 let radius = cell_width / 4.0;
@@ -45,11 +46,11 @@ pub fn elastic_cut_demo(backend: &WebGpu, app_state: &mut AppState) -> PhysicsCo
     };
 
     let params = SimulationParams {
-        gravity: vector![0.0, -9.81, 0.0] * app_state.gravity_factor,
+        gravity: vec3(0.0, -9.81, 0.0) * app_state.gravity_factor,
         dt: 1.0 / 60.0,
     };
 
-    let rb = RigidBodyBuilder::fixed().translation(vector![0.0, -4.0, 0.0]);
+    let rb = RigidBodyBuilder::fixed().translation(vec3(0.0, -4.0, 0.0));
     let rb_handle = rapier_data.bodies.insert(rb);
     let co = ColliderBuilder::cuboid(100.0, 1.0, 100.0);
     rapier_data
@@ -59,12 +60,13 @@ pub fn elastic_cut_demo(backend: &WebGpu, app_state: &mut AppState) -> PhysicsCo
     // TODO: use only two rectangle per cutting tool.
     //       We can’t right now since we don’t really sample the triangles.
     for k in 0..3 {
-        let heights = DMatrix::zeros(10, 10);
-        let heightfield = HeightField::new(heights, vector![35.0, 1.0, 10.0]);
+        let heights = Array2::zeros(10, 10);
+        let heightfield = HeightField::new(heights, vec3(35.0, 1.0, 10.0));
         let (mut vtx, idx) = heightfield.to_trimesh();
+        let rotation = Quat::from_scaled_axis(vec3(1.3, 0.0, 0.0));
+        let translation = vec3(0.0, 10.0, k as f32 * 10.0 - 10.0);
         vtx.iter_mut().for_each(|pt| {
-            *pt = Isometry3::rotation(vector![1.3, 0.0, 0.0]) * *pt
-                + vector![0.0, 10.0, k as f32 * 10.0 - 10.0]
+            *pt = rotation * *pt + translation;
         });
         let rb = RigidBodyBuilder::fixed();
         let rb_handle = rapier_data.bodies.insert(rb);
diff --git a/crates/slosh3d/examples/heightfield3.rs b/crates/slosh3d/examples/heightfield3.rs
index ee39f2e..5ad0346 100644
--- a/crates/slosh3d/examples/heightfield3.rs
+++ b/crates/slosh3d/examples/heightfield3.rs
@@ -1,7 +1,8 @@
 use slosh_testbed3d::{RapierData, slosh};
 
-use nalgebra::{DMatrix, point, vector};
+use glam::vec3;
 use rapier3d::geometry::{HeightField, TriMeshFlags};
+use rapier3d::parry::utils::Array2;
 use rapier3d::prelude::{ColliderBuilder, RigidBodyBuilder};
 use slang_hal::backend::WebGpu;
 use slosh::{
@@ -24,11 +25,11 @@ pub fn heightfield_demo(backend: &WebGpu, app_state: &mut AppState) -> PhysicsCo
     for i in 0..nxz {
         for j in 0..100 {
             for k in 0..nxz {
-                let position = point![
+                let position = vec3(
                     i as f32 + 0.5 - nxz as f32 / 2.0,
                     j as f32 + 0.5 + 14.0,
-                    k as f32 + 0.5 - nxz as f32 / 2.0
-                ] * cell_width
+                    k as f32 + 0.5 - nxz as f32 / 2.0,
+                ) * cell_width
                     / 2.0;
                 let density = 2700.0;
                 let radius = cell_width / 4.0;
@@ -45,14 +46,14 @@ pub fn heightfield_demo(backend: &WebGpu, app_state: &mut AppState) -> PhysicsCo
     };
 
     let params = SimulationParams {
-        gravity: vector![0.0, -9.81, 0.0] * app_state.gravity_factor,
+        gravity: vec3(0.0, -9.81, 0.0) * app_state.gravity_factor,
         dt: 1.0 / 60.0,
     };
 
-    let heights = DMatrix::from_fn(200, 200, |i, j| {
+    let heights = Array2::from_fn(200, 200, |i, j| {
         (i as f32 / 10.0).sin() * (j as f32 / 10.0).cos()
     });
-    let heightfield = HeightField::new(heights, vector![100.0, 5.0, 100.0]);
+    let heightfield = HeightField::new(heights, vec3(100.0, 5.0, 100.0));
     let (vtx, idx) = heightfield.to_trimesh();
     let rb = RigidBodyBuilder::fixed();
     let rb_handle = rapier_data.bodies.insert(rb);
diff --git a/crates/slosh3d/examples/sand3.rs b/crates/slosh3d/examples/sand3.rs
index 817faa6..10dcc68 100644
--- a/crates/slosh3d/examples/sand3.rs
+++ b/crates/slosh3d/examples/sand3.rs
@@ -1,6 +1,6 @@
 use slosh_testbed3d::{RapierData, slosh};
 
-use nalgebra::{point, vector};
+use glam::vec3;
 use rapier3d::prelude::{ColliderBuilder, RigidBodyBuilder};
 use slang_hal::backend::WebGpu;
 use slosh::{
@@ -29,11 +29,11 @@ pub fn sand_demo(backend: &WebGpu, app_state: &mut AppState) -> PhysicsContext {
     for i in 0..nxz {
         for j in 0..100 {
             for k in 0..nxz {
-                let position = point![
+                let position = vec3(
                     i as f32 + 0.5 - nxz as f32 / 2.0,
                     j as f32 + 0.5 + 10.0,
-                    k as f32 + 0.5 - nxz as f32 / 2.0
-                ] * cell_width
+                    k as f32 + 0.5 - nxz as f32 / 2.0,
+                ) * cell_width
                     / 2.0;
                 let radius = cell_width / 4.0;
                 let model = ParticleModel::sand(YOUNG_MODULUS, POISSON_RATIO);
@@ -73,36 +73,36 @@ pub fn sand_demo(backend: &WebGpu, app_state: &mut AppState) -> PhysicsContext {
     };
 
     let params = SimulationParams {
-        gravity: vector![0.0, -9.81, 0.0] * app_state.gravity_factor,
+        gravity: vec3(0.0, -9.81, 0.0) * app_state.gravity_factor,
         dt: 1.0 / 60.0,
     };
 
-    let rb = RigidBodyBuilder::fixed().translation(vector![0.0, -4.0, 0.0]);
+    let rb = RigidBodyBuilder::fixed().translation(vec3(0.0, -4.0, 0.0));
     let rb_handle = rapier_data.bodies.insert(rb);
     let co = ColliderBuilder::cuboid(100.0, 4.0, 100.0);
     rapier_data
         .colliders
         .insert_with_parent(co, rb_handle, &mut rapier_data.bodies);
 
-    let rb = RigidBodyBuilder::fixed().translation(vector![0.0, 5.0, -35.0]);
+    let rb = RigidBodyBuilder::fixed().translation(vec3(0.0, 5.0, -35.0));
     let rb_handle = rapier_data.bodies.insert(rb);
     let co = ColliderBuilder::cuboid(35.0, 5.0, 0.5);
     rapier_data
         .colliders
         .insert_with_parent(co, rb_handle, &mut rapier_data.bodies);
-    let rb = RigidBodyBuilder::fixed().translation(vector![0.0, 5.0, 35.0]);
+    let rb = RigidBodyBuilder::fixed().translation(vec3(0.0, 5.0, 35.0));
     let rb_handle = rapier_data.bodies.insert(rb);
     let co = ColliderBuilder::cuboid(35.0, 5.0, 0.5);
     rapier_data
         .colliders
         .insert_with_parent(co, rb_handle, &mut rapier_data.bodies);
-    let rb = RigidBodyBuilder::fixed().translation(vector![-35.0, 5.0, 0.0]);
+    let rb = RigidBodyBuilder::fixed().translation(vec3(-35.0, 5.0, 0.0));
     let rb_handle = rapier_data.bodies.insert(rb);
     let co = ColliderBuilder::cuboid(0.5, 5.0, 35.0);
     rapier_data
         .colliders
         .insert_with_parent(co, rb_handle, &mut rapier_data.bodies);
-    let rb = RigidBodyBuilder::fixed().translation(vector![35.0, 5.0, 0.0]);
+    let rb = RigidBodyBuilder::fixed().translation(vec3(35.0, 5.0, 0.0));
     let rb_handle = rapier_data.bodies.insert(rb);
     let co = ColliderBuilder::cuboid(0.5, 5.0, 35.0);
     rapier_data
@@ -110,9 +110,9 @@ pub fn sand_demo(backend: &WebGpu, app_state: &mut AppState) -> PhysicsContext {
         .insert_with_parent(co, rb_handle, &mut rapier_data.bodies);
 
     let rb = RigidBodyBuilder::kinematic_velocity_based()
-        .translation(vector![0.0, 2.0, 0.0])
-        .rotation(vector![0.0, 0.0, -0.5])
-        .angvel(vector![0.0, -1.0, 0.0]);
+        .translation(vec3(0.0, 2.0, 0.0))
+        .rotation(vec3(0.0, 0.0, -0.5))
+        .angvel(vec3(0.0, -1.0, 0.0));
     let rb_handle = rapier_data.bodies.insert(rb);
     let co = ColliderBuilder::cuboid(0.5, 2.0, 30.0);
     rapier_data
diff --git a/crates/slosh_testbed2d/Cargo.toml b/crates/slosh_testbed2d/Cargo.toml
index 6a1ceff..13ec6f5 100644
--- a/crates/slosh_testbed2d/Cargo.toml
+++ b/crates/slosh_testbed2d/Cargo.toml
@@ -18,12 +18,13 @@ required-features = ["dim2"]
 dim2 = []
 default = ["dim2"]
 
-comptime = [ "slosh2d/comptime" ]
-runtime = [ "slosh2d/runtime" ]
+comptime = ["slosh2d/comptime"]
+runtime = ["slosh2d/runtime"]
 
-webgpu = [ "slosh2d/webgpu" ]
+webgpu = ["slosh2d/webgpu"]
 
 [dependencies]
+glam = { workspace = true }
 nalgebra = { workspace = true, features = ["rand"] }
 wgpu = { workspace = true }
 bytemuck = { workspace = true }
@@ -41,7 +42,7 @@ slosh2d = { version = "0.6", path = "../slosh2d" }
 regex = "1"
 web-time = "1"
 
-kiss3d = { version = "0.37", features = [ "vertex_index_u32", "egui" ] }
+kiss3d = { version = "0.37", features = ["vertex_index_u32", "egui"] }
 
 
 [build-dependencies]
diff --git a/crates/slosh_testbed3d/Cargo.toml b/crates/slosh_testbed3d/Cargo.toml
index 076bc4a..af50ded 100644
--- a/crates/slosh_testbed3d/Cargo.toml
+++ b/crates/slosh_testbed3d/Cargo.toml
@@ -18,12 +18,13 @@ required-features = ["dim3"]
 dim3 = []
 default = ["dim3"]
 
-comptime = [ "slosh3d/comptime" ]
-runtime = [ "slosh3d/runtime" ]
+comptime = ["slosh3d/comptime"]
+runtime = ["slosh3d/runtime"]
 
-webgpu = [ "slosh3d/webgpu" ]
+webgpu = ["slosh3d/webgpu"]
 
 [dependencies]
+glam = { workspace = true }
 nalgebra = { workspace = true, features = ["rand"] }
 wgpu = { workspace = true }
 bytemuck = { workspace = true }
@@ -41,7 +42,7 @@ slosh3d = { version = "0.6", path = "../slosh3d" }
 regex = "1"
 web-time = "1"
 
-kiss3d = { version = "0.37", features = [ "vertex_index_u32", "egui" ] }
+kiss3d = { version = "0.37", features = ["vertex_index_u32", "egui"] }
 
 
 [build-dependencies]
diff --git a/src/grid/grid.rs b/src/grid/grid.rs
index 78bf651..f4250f5 100644
--- a/src/grid/grid.rs
+++ b/src/grid/grid.rs
@@ -2,7 +2,7 @@
 
 use crate::grid::prefix_sum::{PrefixSumWorkspace, WgPrefixSum};
 use crate::grid::sort::WgSort;
-use crate::math::{Point, Vector};
+use crate::math::Vector;
 use crate::solver::{GpuParticleModelData, GpuParticles, GpuRigidParticles, ParticlePosition};
 use bytemuck::{Pod, Zeroable};
 use encase::ShaderType;
@@ -44,7 +44,7 @@ struct GridArgs<'a, B: Backend> {
     particles_pos: &'a GpuVector<ParticlePosition, B>,
     particles_len: &'a GpuScalar<u32, B>,
     active_blocks: &'a GpuVector<GpuActiveBlockHeader, B>,
-    rigid_particles_pos: &'a GpuVector<Point<f32>, B>,
+    rigid_particles_pos: &'a GpuVector<Vector, B>,
     rigid_particle_needs_block: &'a GpuVector<u32, B>,
     sorted_particle_ids: &'a GpuVector<u32, B>,
     particle_node_linked_lists: &'a GpuVector<u32, B>,
@@ -240,8 +240,8 @@ pub struct GpuGridMetadata {
 #[derive(Copy, Clone, PartialEq, ShaderType)]
 #[repr(C)]
 pub struct GpuGridNode {
-    momentum_velocity_mass: nalgebra::Vector4<f32>,
-    momentum_velocity_mass_incompatible: nalgebra::Vector4<f32>,
+    momentum_velocity_mass: glam::Vec4,
+    momentum_velocity_mass_incompatible: glam::Vec4,
     cdf: GpuGridNodeCdf,
 }
 
@@ -253,9 +253,9 @@ pub struct GpuGridNode {
 #[repr(C)]
 pub struct BlockVirtualId {
     #[cfg(feature = "dim2")]
-    id: nalgebra::Vector2<i32>,
+    id: glam::IVec2,
     #[cfg(feature = "dim3")]
-    id: nalgebra::Vector4<i32>, // Vector3 with padding.
+    id: glam::IVec4, // IVec3 with padding.
 }
 
 /// Hash map entry mapping virtual block IDs to physical storage indices.
@@ -269,13 +269,13 @@ pub struct GpuGridHashMapEntry {
     #[cfg(feature = "dim2")]
     pad0: u32,
     #[cfg(feature = "dim3")]
-    pad0: nalgebra::Vector3<u32>,
+    pad0: glam::UVec3,
     key: BlockVirtualId,
     value: u32,
     #[cfg(feature = "dim2")]
     pad1: u32,
     #[cfg(feature = "dim3")]
-    pad1: nalgebra::Vector3<u32>,
+    pad1: glam::UVec3,
 }
 
 impl Default for GpuGridHashMapEntry {
@@ -331,9 +331,9 @@ pub struct GpuActiveBlockHeader {
 #[derive(Copy, Clone, Default, ShaderType)]
 pub struct GpuNodeCollision {
     /// Outward collision normal at the node.
-    normal: Vector<f32>,
+    normal: Vector,
     /// Closest point on the collider surface.
-    point: Vector<f32>,
+    point: Vector,
     /// ID of the colliding shape.
     collider_id: u32,
     /// Whether a collision was detected (stored as a uint for layout reasons).
diff --git a/src/grid/prefix_sum.rs b/src/grid/prefix_sum.rs
index ae1312a..2a6c544 100644
--- a/src/grid/prefix_sum.rs
+++ b/src/grid/prefix_sum.rs
@@ -1,6 +1,5 @@
 //! Parallel prefix sum (scan) implementation for GPU.
 
-use nalgebra::DVector;
 use slang_hal::backend::Backend;
 use slang_hal::function::GpuFunction;
 use slang_hal::{BufferUsages, Shader, ShaderArgs};
@@ -95,7 +94,7 @@ impl<B: Backend> WgPrefixSum<B> {
     /// CPU implementation of the prefix sum for testing/validation.
     ///
     /// Applies the same algorithm as the GPU version but on CPU.
-    pub fn eval_cpu(&self, v: &mut DVector<u32>) {
+    pub fn eval_cpu(&self, v: &mut [u32]) {
         for i in 0..v.len() - 1 {
             v[i + 1] += v[i];
         }
@@ -158,7 +157,7 @@ impl<B: Backend> PrefixSumWorkspace<B> {
             while stage_len != 1 {
                 let buffer = GpuTensor::vector(
                     backend,
-                    DVector::<u32>::zeros(stage_len as usize),
+                    vec![0u32; stage_len as usize],
                     BufferUsages::STORAGE,
                 )?;
                 self.stages.push(PrefixSumStage {
@@ -172,11 +171,7 @@ impl<B: Backend> PrefixSumWorkspace<B> {
             // The last stage always has only 1 element.
             self.stages.push(PrefixSumStage {
                 capacity: 1,
-                buffer: GpuTensor::vector(
-                    backend,
-                    DVector::<u32>::zeros(1),
-                    BufferUsages::STORAGE,
-                )?,
+                buffer: GpuTensor::vector(backend, vec![0u32; 1], BufferUsages::STORAGE)?,
             });
             self.num_stages = self.stages.len();
         } else if self.stages[0].buffer.len() as u32 != stage_len {
diff --git a/src/grid/sort.rs b/src/grid/sort.rs
index 92020f6..1b3329b 100644
--- a/src/grid/sort.rs
+++ b/src/grid/sort.rs
@@ -1,7 +1,7 @@
 //! Particle sorting kernels for spatial acceleration.
 
 use crate::grid::grid::{GpuGrid, GpuGridHashMapEntry, GpuGridMetadata};
-use crate::math::Point;
+use crate::math::Vector;
 use crate::solver::GpuRigidParticles;
 use slang_hal::backend::Backend;
 use slang_hal::function::GpuFunction;
@@ -39,7 +39,7 @@ struct SortArgs<'a, B: Backend> {
     grid: &'a GpuScalar<GpuGridMetadata, B>,
     hmap_entries: &'a GpuScalar<GpuGridHashMapEntry, B>,
     rigid_nodes_linked_lists: &'a GpuScalar<[u32; 2], B>,
-    rigid_particles_pos: &'a GpuScalar<Point<f32>, B>,
+    rigid_particles_pos: &'a GpuScalar<Vector, B>,
     rigid_particle_node_linked_lists: &'a GpuScalar<u32, B>,
 }
 
diff --git a/src/lib.rs b/src/lib.rs
index b0dec02..ee509d6 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -108,10 +108,10 @@ pub mod math {
 
     /// Angular inertia type for 2D simulations (scalar).
     #[cfg(feature = "dim2")]
-    pub type AngularInertia<N> = N;
+    pub type AngularInertia = f32;
     /// Angular inertia type for 3D simulations (3x3 matrix).
     #[cfg(feature = "dim3")]
-    pub type AngularInertia<N> = rapier::na::Matrix3<N>;
+    pub type AngularInertia = glam::Mat3;
 }
 
 /// Re-exports of commonly used dependencies for convenience.
diff --git a/src/pipeline.rs b/src/pipeline.rs
index 0f38cfa..aaebd47 100644
--- a/src/pipeline.rs
+++ b/src/pipeline.rs
@@ -175,7 +175,8 @@ impl<B: Backend, GpuModel: GpuParticleModelData> MpmPipelineHooks<B, GpuModel> f
 pub struct MpmData<B: Backend, GpuModel: GpuParticleModelData> {
     /// The simulation timestep.
     pub base_dt: f32,
-    pub gravity: Vector<f32>,
+    /// Gravitational acceleration vector (m/s²).
+    pub gravity: Vector,
     /// Global simulation parameters (gravity, timestep).
     pub sim_params: GpuSimulationParams<B>,
     /// Spatial grid for momentum transfer.
diff --git a/src/rbd/dynamics/body.rs b/src/rbd/dynamics/body.rs
index f8687c6..df9be3e 100644
--- a/src/rbd/dynamics/body.rs
+++ b/src/rbd/dynamics/body.rs
@@ -1,9 +1,8 @@
 //! Rigid-body definition and set.
 
-use crate::math::{AngularInertia, GpuSim};
+use crate::math::{AngVector, AngularInertia, GpuSim, Vector};
 use crate::rbd::shapes::{GpuShape, ShapeBuffers};
 use rapier::geometry::ColliderHandle;
-use rapier::math::{AngVector, Point, Vector};
 use rapier::prelude::MassProperties;
 use rapier::{
     dynamics::{RigidBodyHandle, RigidBodySet},
@@ -17,9 +16,9 @@ use stensor::tensor::GpuTensor;
 /// Linear and angular forces with a layout compatible with the corresponding WGSL struct.
 pub struct GpuForce {
     /// The linear part of the force.
-    pub linear: Vector<f32>,
+    pub linear: Vector,
     /// The angular part of the force (aka. the torque).
-    pub angular: AngVector<f32>,
+    pub angular: AngVector,
 }
 
 #[derive(Copy, Clone, PartialEq, Default, encase::ShaderType)]
@@ -27,9 +26,9 @@ pub struct GpuForce {
 /// Linear and angular velocities with a layout compatible with the corresponding WGSL struct.
 pub struct GpuVelocity {
     /// The linear (translational) velocity.
-    pub linear: Vector<f32>,
+    pub linear: Vector,
     /// The angular (rotational) velocity.
-    pub angular: AngVector<f32>,
+    pub angular: AngVector,
 }
 
 #[derive(Copy, Clone, PartialEq, encase::ShaderType)]
@@ -37,11 +36,11 @@ pub struct GpuVelocity {
 /// Rigid-body mass-properties, with a layout compatible with the corresponding WGSL struct.
 pub struct GpuMassProperties {
     /// The inverse angular inertia tensor.
-    pub inv_inertia: AngularInertia<f32>,
+    pub inv_inertia: AngularInertia,
     /// The inverse mass.
-    pub inv_mass: Vector<f32>,
+    pub inv_mass: Vector,
     /// The center-of-mass.
-    pub com: Vector<f32>, // ShaderType isn’t implemented for Point
+    pub com: Vector,
 }
 
 impl From<MassProperties> for GpuMassProperties {
@@ -51,8 +50,8 @@ impl From<MassProperties> for GpuMassProperties {
             inv_inertia: props.inv_principal_inertia,
             #[cfg(feature = "dim3")]
             inv_inertia: props.reconstruct_inverse_inertia_matrix(),
-            inv_mass: Vector::repeat(props.inv_mass),
-            com: props.local_com.coords,
+            inv_mass: Vector::splat(props.inv_mass),
+            com: props.local_com,
         }
     }
 }
@@ -64,9 +63,9 @@ impl Default for GpuMassProperties {
             #[cfg(feature = "dim2")]
             inv_inertia: 1.0,
             #[cfg(feature = "dim3")]
-            inv_inertia: AngularInertia::identity(),
-            inv_mass: Vector::repeat(1.0),
-            com: Vector::zeros(),
+            inv_inertia: AngularInertia::IDENTITY,
+            inv_mass: Vector::splat(1.0),
+            com: Vector::ZERO,
         }
     }
 }
@@ -82,11 +81,11 @@ pub struct GpuBodySet<B: Backend> {
     // TODO: support other shape types.
     // TODO: support a shape with a shift relative to the body.
     pub(crate) shapes: GpuTensor<GpuShape, B>,
-    pub(crate) shapes_local_vertex_buffers: GpuTensor<Point<f32>, B>,
-    pub(crate) shapes_vertex_buffers: GpuTensor<Point<f32>, B>,
+    pub(crate) shapes_local_vertex_buffers: GpuTensor<Vector, B>,
+    pub(crate) shapes_vertex_buffers: GpuTensor<Vector, B>,
     pub(crate) shapes_vertex_collider_id: GpuTensor<u32, B>, // NOTE: this is a bit of a hack for wgsparkl
     /// Vertex buffer for trimesh collision (BVH AABBs, vertices, pseudo-normals).
-    pub(crate) shapes_collision_vertices: GpuTensor<Point<f32>, B>,
+    pub(crate) shapes_collision_vertices: GpuTensor<Vector, B>,
     /// Index buffer for trimesh collision (BVH topology, triangle indices).
     pub(crate) shapes_collision_indices: GpuTensor<u32, B>,
 }
@@ -113,7 +112,7 @@ impl Default for BodyDesc {
             mprops: Default::default(),
             vel: Default::default(),
             pose: Default::default(),
-            shape: GpuShape::cuboid(Vector::repeat(0.5)),
+            shape: GpuShape::cuboid(Vector::splat(0.5)),
         }
     }
 }
@@ -206,14 +205,14 @@ impl<B: Backend> GpuBodySet<B> {
             let two_ways_coupling = rb.is_dynamic() && coupling.mode == BodyCoupling::TwoWays;
             let desc = BodyDesc {
                 vel: GpuVelocity {
-                    linear: *rb.linvel(),
+                    linear: rb.linvel(),
                     #[allow(clippy::clone_on_copy)] // Needed for 2D/3D switch.
                     angular: rb.angvel().clone(),
                 },
                 #[cfg(feature = "dim2")]
-                pose: (*rb.position()).into(),
+                pose: GpuSim::from(rapier::na::Isometry2::from(*rb.position())),
                 #[cfg(feature = "dim3")]
-                pose: GpuSim::from_isometry(*rb.position(), 1.0),
+                pose: GpuSim::from_isometry((*rb.position()).into(), 1.0),
                 shape,
                 local_mprops: if two_ways_coupling {
                     rb.mass_properties().local_mprops.into()
@@ -257,7 +256,7 @@ impl<B: Backend> GpuBodySet<B> {
         // zero-sized buffer bindings. A single dummy element is harmless because
         // the shader only accesses these buffers through index ranges stored in
         // the per-shape data, so the dummy element is never read.
-        let dummy_pt = Point::origin();
+        let dummy_pt = Vector::ZERO;
         let vertices = if shape_buffers.vertices.is_empty() {
             vec![dummy_pt]
         } else {
@@ -371,7 +370,7 @@ impl<B: Backend> GpuBodySet<B> {
     ///
     /// Contains vertices for polylines and trimeshes in world-space coordinates.
     /// Updated when body poses change.
-    pub fn shapes_vertex_buffers(&self) -> &GpuTensor<Point<f32>, B> {
+    pub fn shapes_vertex_buffers(&self) -> &GpuTensor<Vector, B> {
         &self.shapes_vertex_buffers
     }
 
@@ -379,7 +378,7 @@ impl<B: Backend> GpuBodySet<B> {
     ///
     /// Contains vertices for polylines and trimeshes in body-local coordinates.
     /// These are the original untransformed vertices.
-    pub fn shapes_local_vertex_buffers(&self) -> &GpuTensor<Point<f32>, B> {
+    pub fn shapes_local_vertex_buffers(&self) -> &GpuTensor<Vector, B> {
         &self.shapes_local_vertex_buffers
     }
 
@@ -394,7 +393,7 @@ impl<B: Backend> GpuBodySet<B> {
     /// GPU storage buffer containing collision vertices for trimesh shapes.
     ///
     /// Contains BVH AABBs, mesh vertices, and pseudo-normals in body-local coordinates.
-    pub fn shapes_collision_vertices(&self) -> &GpuTensor<Point<f32>, B> {
+    pub fn shapes_collision_vertices(&self) -> &GpuTensor<Vector, B> {
         &self.shapes_collision_vertices
     }
 
diff --git a/src/rbd/shapes.rs b/src/rbd/shapes.rs
index eb91cb9..126239d 100644
--- a/src/rbd/shapes.rs
+++ b/src/rbd/shapes.rs
@@ -4,10 +4,10 @@
 //! optimized for GPU computation. It includes conversion utilities from Rapier/Parry
 //! shapes to GPU-friendly formats with vertex buffers.
 
-use crate::rapier::na::{Vector4, vector};
+use glam::Vec4;
 use rapier::geometry::{Shape, ShapeType, TypedShape};
 
-use crate::math::{Point, Vector};
+use crate::math::Vector;
 
 /// GPU shape type identifiers.
 ///
@@ -43,9 +43,9 @@ pub struct ShapeBuffers {
     /// Vertex positions for all complex shapes.
     ///
     /// Polyline and trimesh shapes store references to ranges within this buffer.
-    pub vertices: Vec<Point<f32>>,
+    pub vertices: Vec<Vector>,
     /// Vertex buffer for trimesh collision (BVH AABBs, vertices, pseudo-normals).
-    pub collision_vertices: Vec<Point<f32>>,
+    pub collision_vertices: Vec<Vector>,
     /// Index buffer for trimesh collision (BVH topology, triangle indices).
     pub collision_indices: Vec<u32>,
 }
@@ -64,8 +64,8 @@ pub struct ShapeBuffers {
 #[derive(Copy, Clone, Debug, bytemuck::Pod, bytemuck::Zeroable)]
 #[repr(C)]
 pub struct GpuShape {
-    a: Vector4<f32>,
-    b: Vector4<f32>,
+    a: Vec4,
+    b: Vec4,
 }
 
 impl GpuShape {
@@ -82,8 +82,8 @@ impl GpuShape {
     pub fn ball(radius: f32) -> Self {
         let tag = f32::from_bits(GpuShapeType::Ball as u32);
         Self {
-            a: vector![radius, 0.0, 0.0, tag],
-            b: vector![0.0, 0.0, 0.0, 0.0],
+            a: Vec4::new(radius, 0.0, 0.0, tag),
+            b: Vec4::ZERO,
         }
     }
 
@@ -98,14 +98,14 @@ impl GpuShape {
     /// # use crate::math::vector;
     /// let cuboid = GpuShape::cuboid(vector![1.0, 2.0, 3.0]);
     /// ```
-    pub fn cuboid(half_extents: Vector<f32>) -> Self {
+    pub fn cuboid(half_extents: Vector) -> Self {
         let tag = f32::from_bits(GpuShapeType::Cuboid as u32);
         Self {
             #[cfg(feature = "dim2")]
-            a: vector![half_extents.x, half_extents.y, 0.0, tag],
+            a: Vec4::new(half_extents.x, half_extents.y, 0.0, tag),
             #[cfg(feature = "dim3")]
-            a: vector![half_extents.x, half_extents.y, half_extents.z, tag],
-            b: vector![0.0, 0.0, 0.0, 0.0],
+            a: Vec4::new(half_extents.x, half_extents.y, half_extents.z, tag),
+            b: Vec4::ZERO,
         }
     }
 
@@ -117,17 +117,17 @@ impl GpuShape {
     /// * `a` - First endpoint of the capsule's central segment
     /// * `b` - Second endpoint of the capsule's central segment
     /// * `radius` - The radius of the capsule's rounded ends
-    pub fn capsule(a: Point<f32>, b: Point<f32>, radius: f32) -> Self {
+    pub fn capsule(a: Vector, b: Vector, radius: f32) -> Self {
         let tag = f32::from_bits(GpuShapeType::Capsule as u32);
         #[cfg(feature = "dim2")]
         return Self {
-            a: vector![a.x, a.y, 0.0, tag],
-            b: vector![b.x, b.y, 0.0, radius],
+            a: Vec4::new(a.x, a.y, 0.0, tag),
+            b: Vec4::new(b.x, b.y, 0.0, radius),
         };
         #[cfg(feature = "dim3")]
         return Self {
-            a: vector![a.x, a.y, a.z, tag],
-            b: vector![b.x, b.y, b.z, radius],
+            a: Vec4::new(a.x, a.y, a.z, tag),
+            b: Vec4::new(b.x, b.y, b.z, radius),
         };
     }
 
@@ -142,8 +142,8 @@ impl GpuShape {
         let rng0 = f32::from_bits(vertex_range[0]);
         let rng1 = f32::from_bits(vertex_range[1]);
         Self {
-            a: vector![rng0, rng1, 0.0, tag],
-            b: vector![0.0, 0.0, 0.0, 0.0],
+            a: Vec4::new(rng0, rng1, 0.0, tag),
+            b: Vec4::ZERO,
         }
     }
 
@@ -157,18 +157,18 @@ impl GpuShape {
     pub fn trimesh(trimesh: &crate::trimesh::GpuTriMesh, vertex_base_id: u32) -> Self {
         let tag = f32::from_bits(GpuShapeType::TriMesh as u32);
         Self {
-            a: vector![
+            a: Vec4::new(
                 f32::from_bits(trimesh.bvh_vtx_root_id),
                 f32::from_bits(trimesh.bvh_idx_root_id),
                 f32::from_bits(trimesh.bvh_node_len),
-                tag
-            ],
-            b: vector![
+                tag,
+            ),
+            b: Vec4::new(
                 f32::from_bits(trimesh.num_triangles),
                 f32::from_bits(trimesh.num_vertices),
                 f32::from_bits(vertex_base_id),
-                0.0
-            ],
+                0.0,
+            ),
         }
     }
 
@@ -181,8 +181,8 @@ impl GpuShape {
     pub fn cone(half_height: f32, radius: f32) -> Self {
         let tag = f32::from_bits(GpuShapeType::Cone as u32);
         Self {
-            a: vector![half_height, radius, 0.0, tag],
-            b: vector![0.0, 0.0, 0.0, 0.0],
+            a: Vec4::new(half_height, radius, 0.0, tag),
+            b: Vec4::ZERO,
         }
     }
 
@@ -195,8 +195,8 @@ impl GpuShape {
     pub fn cylinder(half_height: f32, radius: f32) -> Self {
         let tag = f32::from_bits(GpuShapeType::Cylinder as u32);
         Self {
-            a: vector![half_height, radius, 0.0, tag],
-            b: vector![0.0, 0.0, 0.0, 0.0],
+            a: Vec4::new(half_height, radius, 0.0, tag),
+            b: Vec4::ZERO,
         }
     }
 
diff --git a/src/sampling/sample_polyline.rs b/src/sampling/sample_polyline.rs
index 71ff4fe..45248d2 100644
--- a/src/sampling/sample_polyline.rs
+++ b/src/sampling/sample_polyline.rs
@@ -1,11 +1,12 @@
+use crate::math::Vector;
 use encase::ShaderType;
-use nalgebra::{Point2, Vector2, vector};
+use glam::UVec2;
 use rapier::geometry::{Polyline, Segment};
 
 #[derive(Copy, Clone, Debug, ShaderType)]
 #[repr(C)]
 pub struct GpuSampleIds {
-    pub segment: Vector2<u32>,
+    pub segment: UVec2,
     pub collider: u32,
 }
 
@@ -19,7 +20,7 @@ pub struct SamplingParams {
 
 #[derive(Default, Clone)]
 pub struct SamplingBuffers {
-    pub samples: Vec<Point2<f32>>,
+    pub samples: Vec<Vector>,
     pub samples_ids: Vec<GpuSampleIds>,
 }
 
@@ -34,7 +35,7 @@ pub fn sample_polyline(
             polyline.vertices()[seg_idx[1] as usize],
         );
         let sample_id = GpuSampleIds {
-            segment: vector![params.base_vid + seg_idx[0], params.base_vid + seg_idx[1]],
+            segment: UVec2::new(params.base_vid + seg_idx[0], params.base_vid + seg_idx[1]),
             collider: params.collider_id,
         };
         buffers.samples.push(seg.a);
@@ -47,7 +48,7 @@ pub fn sample_polyline(
                     break;
                 }
 
-                buffers.samples.push(seg.a + *dir * shift);
+                buffers.samples.push(seg.a + dir * shift);
                 buffers.samples_ids.push(sample_id);
             }
 
diff --git a/src/sampling/sample_trimesh.rs b/src/sampling/sample_trimesh.rs
index ad06e72..165d493 100644
--- a/src/sampling/sample_trimesh.rs
+++ b/src/sampling/sample_trimesh.rs
@@ -1,5 +1,6 @@
+use crate::math::Vector;
 use encase::ShaderType;
-use nalgebra::{Point3, Vector3, vector};
+use glam::UVec3;
 use rapier::geometry::{Segment, TriMesh, Triangle};
 use std::collections::HashSet;
 
@@ -9,13 +10,13 @@ const EPS: f32 = 1.0e-5;
 
 pub struct TriangleSample {
     pub triangle_id: u32,
-    pub point: Point3<f32>,
+    pub point: Vector,
 }
 
 #[derive(Copy, Clone, Debug, ShaderType)]
 #[repr(C)]
 pub struct GpuSampleIds {
-    pub triangle: Vector3<u32>,
+    pub triangle: UVec3,
     pub collider: u32,
 }
 
@@ -28,8 +29,8 @@ pub struct SamplingParams {
 
 #[derive(Default, Clone)]
 pub struct SamplingBuffers {
-    pub local_samples: Vec<Point3<f32>>,
-    pub samples: Vec<Point3<f32>>,
+    pub local_samples: Vec<Vector>,
+    pub samples: Vec<Vector>,
     pub samples_ids: Vec<GpuSampleIds>,
 }
 
@@ -39,11 +40,11 @@ pub fn sample_trimesh(trimesh: &TriMesh, params: &SamplingParams, buffers: &mut
     for sample in samples {
         let tri_idx = trimesh.indices()[sample.triangle_id as usize];
         let sample_id = GpuSampleIds {
-            triangle: vector![
+            triangle: UVec3::new(
                 params.base_vid + tri_idx[0],
                 params.base_vid + tri_idx[1],
-                params.base_vid + tri_idx[2]
-            ],
+                params.base_vid + tri_idx[2],
+            ),
             collider: params.collider_id,
         };
         buffers.local_samples.push(sample.point);
@@ -61,7 +62,7 @@ pub fn sample_trimesh(trimesh: &TriMesh, params: &SamplingParams, buffers: &mut
 /// Samples a triangle mesh with a set of points such that at least one point is generated
 /// inside each cell on a grid on the x-y plane with cells sized by `xy_spacing`.
 pub fn sample_mesh(
-    vertices: &[Point3<f32>],
+    vertices: &[Vector],
     indices: &[[u32; 3]],
     xy_spacing: f32,
 ) -> Vec<TriangleSample> {
@@ -117,7 +118,7 @@ pub fn sample_edge(
     triangle_id: u32,
 ) {
     let ab = edge.b - edge.a;
-    let edge_length = ab.norm();
+    let edge_length = ab.length();
 
     if edge_length > EPS {
         let edge_dir = ab / edge_length;
@@ -153,9 +154,9 @@ pub fn sample_triangle(
     triangle_id: u32,
 ) {
     // select the longest edge as the base
-    let distance_ab = nalgebra::distance(&triangle.b, &triangle.a);
-    let distance_bc = nalgebra::distance(&triangle.c, &triangle.b);
-    let distance_ca = nalgebra::distance(&triangle.a, &triangle.c);
+    let distance_ab = triangle.b.distance(triangle.a);
+    let distance_bc = triangle.c.distance(triangle.b);
+    let distance_ca = triangle.a.distance(triangle.c);
     let max = distance_ab.max(distance_bc).max(distance_ca);
 
     let triangle = if max == distance_bc {
@@ -176,7 +177,7 @@ pub fn sample_triangle(
 
     let ac = triangle.c - triangle.a;
     let base = triangle.b - triangle.a;
-    let base_length = base.norm();
+    let base_length = base.length();
     let base_dir = base / base_length;
 
     // Adjust the spacing so it matches the required spacing on the x-y plane.
@@ -194,7 +195,7 @@ pub fn sample_triangle(
     let base_step = base_dir * spacing;
 
     // Project C on the base AB.
-    let ac_offset_length = ac.dot(&base_dir);
+    let ac_offset_length = ac.dot(base_dir);
     let bc_offset_length = base_length - ac_offset_length;
 
     if ac_offset_length < EPS || bc_offset_length < EPS || base_length < EPS {
@@ -203,7 +204,7 @@ pub fn sample_triangle(
 
     // Compute the triangle’s height vector.
     let height = ac - base_dir * ac_offset_length;
-    let height_length = height.norm();
+    let height_length = height.length();
     let height_dir = height / height_length;
     // Calculate the tangents.
     let tan_alpha = height_length / ac_offset_length;
@@ -216,8 +217,8 @@ pub fn sample_triangle(
 
         // Compute the height at the current base_position. The point at the
         // end of that height is either in the line (AC) or (BC), whichever is closer.
-        let height_ac = tan_alpha * nalgebra::distance(&triangle.a, &base_position);
-        let height_bc = tan_beta * nalgebra::distance(&triangle.b, &base_position);
+        let height_ac = tan_alpha * triangle.a.distance(base_position);
+        let height_bc = tan_beta * triangle.b.distance(base_position);
         let height_length = height_ac.min(height_bc);
 
         // Calculate the step increment on the height.
@@ -228,7 +229,7 @@ pub fn sample_triangle(
         for j in 1..height_step_count as u32 {
             let particle_position = base_position + (j as f32) * height_step;
 
-            if particle_position.iter().any(|e| !e.is_finite()) {
+            if !particle_position.is_finite() {
                 continue;
             }
 
diff --git a/src/solver/grid_update.rs b/src/solver/grid_update.rs
index ec4095f..0b0a992 100644
--- a/src/solver/grid_update.rs
+++ b/src/solver/grid_update.rs
@@ -10,7 +10,7 @@ use crate::grid::grid::{
     GpuActiveBlockHeader, GpuGrid, GpuGridHashMapEntry, GpuGridMetadata, GpuGridNode,
     GpuNodeCollision,
 };
-use crate::math::{GpuSim, Point};
+use crate::math::{GpuSim, Vector};
 use crate::rbd::dynamics::GpuBodySet;
 use crate::rbd::shapes::GpuShape;
 use crate::solver::params::GpuSimulationParams;
@@ -42,7 +42,7 @@ struct GridUpdateArgs<'a, B: Backend> {
     body_materials: &'a GpuVector<GpuBoundaryCondition, B>,
     collision_shapes: &'a GpuTensor<GpuShape, B>,
     collision_shape_poses: &'a GpuTensor<GpuSim, B>,
-    collision_shape_vtx: &'a GpuTensor<Point<f32>, B>,
+    collision_shape_vtx: &'a GpuTensor<Vector, B>,
     collision_shape_idx: &'a GpuTensor<u32, B>,
     prev_node_collisions: &'a GpuVector<GpuNodeCollision, B>,
     node_collisions: &'a GpuVector<GpuNodeCollision, B>,
diff --git a/src/solver/p2g_cdf.rs b/src/solver/p2g_cdf.rs
index 06cd30b..3d8f48e 100644
--- a/src/solver/p2g_cdf.rs
+++ b/src/solver/p2g_cdf.rs
@@ -3,7 +3,7 @@
 use crate::grid::grid::{
     GpuActiveBlockHeader, GpuGrid, GpuGridHashMapEntry, GpuGridMetadata, GpuGridNode,
 };
-use crate::math::Point;
+use crate::math::Vector;
 use crate::rbd::dynamics::GpuBodySet;
 use crate::sampling::GpuSampleIds;
 use crate::solver::GpuRigidParticles;
@@ -30,7 +30,7 @@ struct P2GCdfArgs<'a, B: Backend> {
     active_blocks: &'a GpuVector<GpuActiveBlockHeader, B>,
     rigid_nodes_linked_lists: &'a GpuVector<[u32; 2], B>,
     particle_node_linked_lists: &'a GpuVector<u32, B>,
-    collider_vertices: &'a GpuVector<Point<f32>, B>,
+    collider_vertices: &'a GpuVector<Vector, B>,
     rigid_particle_indices: &'a GpuVector<GpuSampleIds, B>,
     nodes: &'a GpuVector<GpuGridNode, B>,
 }
diff --git a/src/solver/params.rs b/src/solver/params.rs
index ab0e51e..59bc317 100644
--- a/src/solver/params.rs
+++ b/src/solver/params.rs
@@ -1,3 +1,4 @@
+use crate::math::Vector;
 use bytemuck::{Pod, Zeroable};
 use slang_hal::{BufferUsages, backend::Backend};
 use stensor::tensor::{GpuScalar, GpuTensor};
@@ -11,14 +12,10 @@ use stensor::tensor::{GpuScalar, GpuTensor};
 #[repr(C)]
 pub struct SimulationParams {
     /// Gravitational acceleration vector (m/s²).
-    #[cfg(feature = "dim2")]
-    pub gravity: nalgebra::Vector2<f32>,
+    pub gravity: Vector,
     /// Padding for GPU alignment (2D only).
     #[cfg(feature = "dim2")]
     pub padding: f32,
-    /// Gravitational acceleration vector (m/s²).
-    #[cfg(feature = "dim3")]
-    pub gravity: nalgebra::Vector3<f32>,
     /// Simulation timestep duration (seconds).
     pub dt: f32,
 }
diff --git a/src/solver/particle.rs b/src/solver/particle.rs
index 81d8860..952be49 100644
--- a/src/solver/particle.rs
+++ b/src/solver/particle.rs
@@ -1,4 +1,4 @@
-use crate::math::{Matrix, Point, Vector};
+use crate::math::{Matrix, Vector};
 use crate::rbd::dynamics::GpuBodySet;
 use crate::rbd::dynamics::body::BodyCouplingEntry;
 use crate::rbd::shapes::ShapeBuffers;
@@ -22,15 +22,15 @@ use stensor::tensor::{GpuScalar, GpuTensor};
 #[repr(C)]
 pub struct ParticleDynamics {
     /// Current velocity (m/s).
-    pub velocity: Vector<f32>,
+    pub velocity: Vector,
     /// Deformation gradient tracking how the particle has deformed from its initial state.
-    pub def_grad: Matrix<f32>,
+    pub def_grad: Matrix,
     /// APIC affine velocity matrix for improved momentum conservation.
-    pub affine: Matrix<f32>,
+    pub affine: Matrix,
     /// Additional force applied indirectly to the particle.
     /// Resets automatically at the `particle_update` stage of the
     /// MPM pipeline.
-    pub force_dt: Vector<f32>,
+    pub force_dt: Vector,
     /// Determinant of velocity gradient (for volume change tracking).
     pub vel_grad_det: f32,
     /// Collision detection field data for rigid body coupling.
@@ -67,10 +67,10 @@ impl ParticleDynamics {
         let exponent = if cfg!(feature = "dim2") { 2 } else { 3 };
         let init_volume = (radius * 2.0).powi(exponent); // NOTE: the particles are square-ish.
         Self {
-            velocity: Vector::zeros(),
-            def_grad: Matrix::identity(),
-            affine: Matrix::zeros(),
-            force_dt: Vector::zeros(),
+            velocity: Vector::ZERO,
+            def_grad: Matrix::IDENTITY,
+            affine: Matrix::ZERO,
+            force_dt: Vector::ZERO,
             vel_grad_det: 0.0,
             init_volume,
             init_radius: radius,
@@ -133,11 +133,11 @@ impl ParticleDynamics {
 #[repr(C)]
 pub struct Kinematics {
     /// APIC affine velocity matrix / velocity gradient.
-    pub affine: Matrix<f32>,
+    pub affine: Matrix,
     /// Current velocity (m/s).
-    pub velocity: Vector<f32>,
+    pub velocity: Vector,
     /// Additional force * dt applied to the particle.
-    pub force_dt: Vector<f32>,
+    pub force_dt: Vector,
     /// Determinant of velocity gradient (for volume change tracking).
     pub vel_grad_det: f32,
     /// Particle mass (kg).
@@ -205,9 +205,9 @@ impl ParticlePhase {
 #[repr(C)]
 pub struct Cdf {
     /// Surface normal of the closest rigid body surface.
-    pub normal: Vector<f32>,
+    pub normal: Vector,
     /// Velocity of the rigid body at the contact point.
-    pub rigid_vel: Vector<f32>,
+    pub rigid_vel: Vector,
     /// Signed distance to the nearest rigid body surface (negative = penetration).
     pub signed_distance: f32,
     /// Index of the rigid body this particle is coupled with.
@@ -225,7 +225,7 @@ pub struct Cdf {
 #[derive(Copy, Clone, Debug)]
 pub struct Particle<Model> {
     /// Spatial position (m).
-    pub position: Point<f32>,
+    pub position: Vector,
     /// Physical state (velocity, deformation, mass, etc.).
     pub dynamics: ParticleDynamics,
     /// Material model defining constitutive behavior.
@@ -241,7 +241,7 @@ impl<Model> Particle<Model> {
     /// * `radius` - Particle radius (used to compute mass and volume)
     /// * `density` - Material density
     /// * `model` - Material model instance
-    pub fn new(position: Point<f32>, radius: f32, density: f32, model: Model) -> Self {
+    pub fn new(position: Vector, radius: f32, density: f32, model: Model) -> Self {
         Particle {
             position,
             dynamics: ParticleDynamics::new(radius, density),
@@ -257,9 +257,9 @@ impl<Model> Particle<Model> {
 /// interact with MPM particles through the CDF (Collision Detection Field).
 pub struct GpuRigidParticles<B: Backend> {
     /// Sample points in local (body-relative) coordinates.
-    pub local_sample_points: GpuTensor<Point<f32>, B>,
+    pub local_sample_points: GpuTensor<Vector, B>,
     /// Sample points transformed to world coordinates.
-    pub sample_points: GpuTensor<Point<f32>, B>,
+    pub sample_points: GpuTensor<Vector, B>,
     /// Bitmask indicating which rigid particles need grid cell blocking.
     pub rigid_particle_needs_block: GpuTensor<u32, B>,
     /// Linked list for spatially sorting rigid particles into grid cells.
@@ -381,18 +381,18 @@ impl<B: Backend> GpuRigidParticles<B> {
     }
 }
 
-/// Particle position type (2D: Point2, 3D: Point4 for alignment).
+/// Particle position type (2D: Vec2, 3D: Vec4 for alignment).
 #[cfg(feature = "dim2")]
-pub type ParticlePosition = Point<f32>;
-/// Particle position type (2D: Point2, 3D: Point4 for alignment).
+pub type ParticlePosition = glam::Vec2;
+/// Particle position type (2D: Vec2, 3D: Vec4 for alignment).
 #[cfg(feature = "dim3")]
-pub type ParticlePosition = nalgebra::Point4<f32>;
+pub type ParticlePosition = glam::Vec4;
 
 struct SoAParticles<GpuModel: GpuParticleModelData> {
     positions: Vec<ParticlePosition>,
     kinematics: Vec<Kinematics>,
     cdf: Vec<Cdf>,
-    def_grad: Vec<Matrix<f32>>,
+    def_grad: Vec<Matrix>,
     properties: Vec<ParticleProperties>,
     models: Vec<GpuModel>,
 }
@@ -402,10 +402,7 @@ impl<GpuModel: GpuParticleModelData> SoAParticles<GpuModel> {
         #[cfg(feature = "dim2")]
         let positions: Vec<_> = particles.iter().map(|p| p.position).collect();
         #[cfg(feature = "dim3")]
-        let positions: Vec<_> = particles
-            .iter()
-            .map(|p| p.position.coords.push(0.0).into())
-            .collect();
+        let positions: Vec<_> = particles.iter().map(|p| p.position.extend(0.0)).collect();
         let kinematics: Vec<_> = particles
             .iter()
             .map(|p| p.dynamics.to_kinematics())
@@ -446,7 +443,7 @@ pub struct GpuParticles<B: Backend, GpuModel: GpuParticleModelData> {
     /// Collision detection field data for rigid body coupling.
     pub cdf: GpuTensor<Cdf, B>,
     /// Deformation gradient matrices.
-    pub def_grad: GpuTensor<Matrix<f32>, B>,
+    pub def_grad: GpuTensor<Matrix, B>,
     /// Static per-particle properties (read-only on GPU).
     pub properties: GpuTensor<ParticleProperties, B>,
     models: GpuTensor<GpuModel, B>,
@@ -487,7 +484,7 @@ impl<B: Backend, GpuModel: GpuParticleModelData> GpuParticles<B, GpuModel> {
                 particles.len() as u32,
                 BufferUsages::STORAGE | BufferUsages::UNIFORM | BufferUsages::COPY_DST,
             )?,
-            positions: GpuTensor::vector(backend, &data.positions, resizeable)?,
+            positions: GpuTensor::vector_encased(backend, &data.positions, resizeable)?,
             kinematics: GpuTensor::vector_encased(backend, &data.kinematics, resizeable)?,
             cdf: GpuTensor::vector_encased(backend, &data.cdf, resizeable)?,
             def_grad: GpuTensor::vector_encased(backend, &data.def_grad, resizeable)?,
diff --git a/src/solver/particle_update.rs b/src/solver/particle_update.rs
index 84b9ad7..216ca52 100644
--- a/src/solver/particle_update.rs
+++ b/src/solver/particle_update.rs
@@ -39,7 +39,7 @@ struct ParticleUpdateArgs<'a, B: Backend, GpuModel: GpuParticleModelData> {
     particles_pos: &'a GpuTensor<ParticlePosition, B>,
     particles_kin: &'a GpuTensor<Kinematics, B>,
     particles_cdf: &'a GpuTensor<Cdf, B>,
-    particles_def_grad: &'a GpuTensor<Matrix<f32>, B>,
+    particles_def_grad: &'a GpuTensor<Matrix, B>,
     particles_props: &'a GpuTensor<ParticleProperties, B>,
     particles_len: &'a GpuScalar<u32, B>,
 }
diff --git a/src/solver/rigid_impulses.rs b/src/solver/rigid_impulses.rs
index ef862e6..c0a846c 100644
--- a/src/solver/rigid_impulses.rs
+++ b/src/solver/rigid_impulses.rs
@@ -1,12 +1,11 @@
 //! Impulse accumulation and application for MPM-rigid body coupling.
 
 use crate::grid::grid::{GpuGrid, GpuGridMetadata};
-use crate::math::GpuSim;
+use crate::math::{AngVector, GpuSim, Vector};
 use crate::rbd::dynamics::{GpuBodySet, GpuMassProperties, GpuVelocity};
 use crate::solver::SimulationParams;
 use crate::solver::params::GpuSimulationParams;
 use encase::ShaderType;
-use rapier::math::{AngVector, Point, Vector};
 use slang_hal::backend::Backend;
 use slang_hal::function::GpuFunction;
 use slang_hal::{BufferUsages, Shader, ShaderArgs};
@@ -32,11 +31,11 @@ pub struct WgRigidImpulses<B: Backend> {
 #[repr(C)]
 pub struct RigidImpulse {
     /// Center of mass (for convenience).
-    pub com: Point<f32>,
+    pub com: Vector,
     /// Linear impulse vector.
-    pub linear: Vector<f32>,
+    pub linear: Vector,
     /// Angular impulse (torque).
-    pub angular: AngVector<f32>,
+    pub angular: AngVector,
 }
 
 /// GPU buffers for storing impulses from MPM to rigid bodies.
diff --git a/src/solver/rigid_particle_update.rs b/src/solver/rigid_particle_update.rs
index ca3d72f..8899ed9 100644
--- a/src/solver/rigid_particle_update.rs
+++ b/src/solver/rigid_particle_update.rs
@@ -1,6 +1,6 @@
 //! Rigid body particle transformation kernels.
 
-use crate::math::{GpuSim, Point};
+use crate::math::{GpuSim, Vector};
 use crate::rbd::dynamics::GpuBodySet;
 use crate::sampling::GpuSampleIds;
 use crate::solver::GpuRigidParticles;
@@ -27,8 +27,8 @@ struct RigidParticleUpdateArgs<'a, B: Backend> {
     vertex_collider_ids: Option<&'a GpuTensor<u32, B>>,
     rigid_particle_indices: Option<&'a GpuTensor<GpuSampleIds, B>>,
     poses: &'a GpuTensor<GpuSim, B>,
-    local_pts: &'a GpuTensor<Point<f32>, B>,
-    world_pts: &'a GpuTensor<Point<f32>, B>,
+    local_pts: &'a GpuTensor<Vector, B>,
+    world_pts: &'a GpuTensor<Vector, B>,
 }
 
 impl<B: Backend> WgRigidParticleUpdate<B> {
diff --git a/src/solver/timestep_bound.rs b/src/solver/timestep_bound.rs
index db6c5c2..8303503 100644
--- a/src/solver/timestep_bound.rs
+++ b/src/solver/timestep_bound.rs
@@ -47,7 +47,7 @@ struct TimestepBoundsArgs<'a, B: Backend, GpuModel: GpuParticleModelData> {
     grid: &'a GpuTensor<GpuGridMetadata, B>,
     particles_model: &'a GpuTensor<GpuModel, B>,
     particles_kin: &'a GpuTensor<Kinematics, B>,
-    particles_def_grad: &'a GpuTensor<Matrix<f32>, B>,
+    particles_def_grad: &'a GpuTensor<Matrix, B>,
     particles_props: &'a GpuTensor<ParticleProperties, B>,
     particles_len: &'a GpuScalar<u32, B>,
     result: &'a GpuScalar<GpuTimestepBounds, B>,
diff --git a/src/trimesh.rs b/src/trimesh.rs
index 118e40b..a6ffc2b 100644
--- a/src/trimesh.rs
+++ b/src/trimesh.rs
@@ -1,8 +1,9 @@
 // TODO: move this to rbd?
 
+use crate::math::Vector;
 use encase::ShaderType;
 use rapier::geometry::TriMesh;
-use rapier::prelude::{DIM, Point};
+use rapier::prelude::DIM;
 
 #[derive(Copy, Clone, ShaderType)]
 pub struct GpuTriMesh {
@@ -20,7 +21,7 @@ pub struct GpuTriMesh {
 
 pub fn convert_trimesh_to_gpu(
     shape: &TriMesh,
-    vertices: &mut Vec<Point<f32>>,
+    vertices: &mut Vec<Vector>,
     indices: &mut Vec<u32>,
 ) -> GpuTriMesh {
     let bvh_vtx_root_id = vertices.len();
@@ -56,7 +57,7 @@ pub fn convert_trimesh_to_gpu(
         .map(|tri| {
             let aabb = tri.local_aabb();
             BvhObject {
-                aabb: bvh::aabb::Aabb::with_bounds(aabb.mins, aabb.maxs),
+                aabb: bvh::aabb::Aabb::with_bounds(aabb.mins.into(), aabb.maxs.into()),
                 node_index: 0,
             }
         })
@@ -64,7 +65,11 @@ pub fn convert_trimesh_to_gpu(
 
     let bvh = bvh::bvh::Bvh::build(&mut objects);
     let flat_bvh = bvh.flatten();
-    vertices.extend(flat_bvh.iter().flat_map(|n| [n.aabb.min, n.aabb.max]));
+    vertices.extend(
+        flat_bvh
+            .iter()
+            .flat_map(|n| [Vector::from(n.aabb.min), Vector::from(n.aabb.max)]),
+    );
     let bvh_node_len = flat_bvh.len();
     indices.extend(
         flat_bvh
@@ -79,12 +84,12 @@ pub fn convert_trimesh_to_gpu(
             .pseudo_normals()
             .expect("trimeshes without pseudo-normals are not supported");
         vertices.extend_from_slice(shape.vertices());
-        vertices.extend(pn.vertices_pseudo_normal.iter().map(|n| Point::from(*n)));
+        vertices.extend_from_slice(&pn.vertices_pseudo_normal);
         assert_eq!(shape.vertices().len(), pn.vertices_pseudo_normal.len());
         vertices.extend(
             pn.edges_pseudo_normal
                 .iter()
-                .flat_map(|n| n.map(Point::from)),
+                .flat_map(|n| n.map(Vector::from)),
         );
     }
     indices.extend(shape.indices().iter().flat_map(|tri| tri.iter().copied()));
diff --git a/src_testbed/data.rs b/src_testbed/data.rs
index 3585b56..b0332fe 100644
--- a/src_testbed/data.rs
+++ b/src_testbed/data.rs
@@ -1,4 +1,5 @@
 use crate::step::SimulationStepResult;
+use glam::Vec4;
 use slang_hal::backend::WebGpu;
 use slosh::pipeline::{MpmData, MpmPipeline};
 use slosh::rapier::prelude::{
@@ -26,7 +27,7 @@ pub struct AppState<GpuModel: GpuParticleModelData = GpuParticleModel> {
     /// the particle count, these override the default index-cycled palette
     /// in the readback color buffer. Scene builders set this field on
     /// `&mut AppState` during `init()` to color particles by material.
-    pub particle_colors: Option<Vec<nalgebra::Vector4<f32>>>,
+    pub particle_colors: Option<Vec<Vec4>>,
     /// Enable the axis-aligned cutting box (3D only). When false, no filtering
     /// is applied and all particles render.
     #[cfg(feature = "dim3")]
@@ -34,16 +35,16 @@ pub struct AppState<GpuModel: GpuParticleModelData = GpuParticleModel> {
     /// Current cutting-box bounds (3D only). Particles outside this box are
     /// culled before being sent to the GPU.
     #[cfg(feature = "dim3")]
-    pub render_aabb_min: nalgebra::Vector3<f32>,
+    pub render_aabb_min: glam::Vec3,
     #[cfg(feature = "dim3")]
-    pub render_aabb_max: nalgebra::Vector3<f32>,
+    pub render_aabb_max: glam::Vec3,
     /// Slider range for the cutting-box sliders (3D only). Scene builders
     /// set this from their particle positions at init time so the sliders
     /// span the actual scene extent.
     #[cfg(feature = "dim3")]
-    pub render_aabb_slider_min: nalgebra::Vector3<f32>,
+    pub render_aabb_slider_min: glam::Vec3,
     #[cfg(feature = "dim3")]
-    pub render_aabb_slider_max: nalgebra::Vector3<f32>,
+    pub render_aabb_slider_max: glam::Vec3,
     /// Optional initial camera eye position (3D only). Set by the scene
     /// builder during `init()` to override the testbed's default arc-ball
     /// camera placement.
diff --git a/src_testbed/lib.rs b/src_testbed/lib.rs
index 67c0676..6a16a89 100644
--- a/src_testbed/lib.rs
+++ b/src_testbed/lib.rs
@@ -28,13 +28,12 @@ use crate::prep_readback::{GpuReadbackData, PrepReadback, RenderMode};
 use kiss3d::egui;
 use kiss3d::planar_camera::Sidescroll;
 use kiss3d::prelude::*;
-#[cfg(feature = "dim3")]
-use nalgebra::Vector3;
 use regex::Regex;
 use slang_hal::backend::{Backend, WebGpu};
 use slang_hal::re_exports::include_dir;
 use slang_hal::BufferUsages;
 use slang_hal::SlangCompiler;
+use slosh::math::Vector;
 use slosh::pipeline::{MpmPipeline, MpmPipelineHooks};
 use slosh::rapier::geometry::Shape;
 use slosh::rapier::geometry::ShapeType;
@@ -125,13 +124,13 @@ impl<GpuModel: GpuParticleModelData> Stage<GpuModel> {
             #[cfg(feature = "dim3")]
             render_aabb_enabled: false,
             #[cfg(feature = "dim3")]
-            render_aabb_min: nalgebra::Vector3::repeat(-100.0),
+            render_aabb_min: Vector::splat(-100.0),
             #[cfg(feature = "dim3")]
-            render_aabb_max: nalgebra::Vector3::repeat(100.0),
+            render_aabb_max: Vector::splat(100.0),
             #[cfg(feature = "dim3")]
-            render_aabb_slider_min: nalgebra::Vector3::repeat(-100.0),
+            render_aabb_slider_min: Vector::splat(-100.0),
             #[cfg(feature = "dim3")]
-            render_aabb_slider_max: nalgebra::Vector3::repeat(100.0),
+            render_aabb_slider_max: Vector::splat(100.0),
             #[cfg(feature = "dim3")]
             initial_camera_eye: None,
             #[cfg(feature = "dim3")]
@@ -247,8 +246,8 @@ impl<GpuModel: GpuParticleModelData> Stage<GpuModel> {
         // inward stay where they are.
         #[cfg(feature = "dim3")]
         if !self.step_result.instances.is_empty() {
-            let mut new_min = nalgebra::Vector3::repeat(f32::INFINITY);
-            let mut new_max = nalgebra::Vector3::repeat(f32::NEG_INFINITY);
+            let mut new_min = Vector::splat(f32::INFINITY);
+            let mut new_max = Vector::splat(f32::NEG_INFINITY);
             for d in &self.step_result.instances {
                 new_min.x = new_min.x.min(d.position.x);
                 new_min.y = new_min.y.min(d.position.y);
@@ -281,7 +280,7 @@ impl<GpuModel: GpuParticleModelData> Stage<GpuModel> {
                 .iter()
                 .map(|d| PlanarInstanceData {
                     position: kiss3d::nalgebra::Point2::new(d.position.x, d.position.y),
-                    color: d.color.into(),
+                    color: d.color.to_array(),
                     #[rustfmt::skip]
                     deformation: kiss3d::nalgebra::Matrix2::new(
                         d.deformation.m11, d.deformation.m12,
@@ -313,7 +312,7 @@ impl<GpuModel: GpuParticleModelData> Stage<GpuModel> {
                             d.position.y,
                             d.position.z,
                         ),
-                        color: d.color.into(),
+                        color: d.color.to_array(),
                         #[rustfmt::skip]
                         deformation: kiss3d::nalgebra::Matrix3::new(
                             d.deformation.m11, d.deformation.m12, d.deformation.m13,
@@ -342,7 +341,7 @@ pub async fn run_with_compiler<GpuModel: GpuParticleModelData>(
     }
     #[cfg(feature = "dim3")]
     {
-        run_with_hooks(compiler, |_, _| Box::new(()), scene_builders, Vector3::y()).await;
+        run_with_hooks(compiler, |_, _| Box::new(()), scene_builders, Vector::Y).await;
     }
 }
 
@@ -350,7 +349,7 @@ pub async fn run_with_hooks<GpuModel: GpuParticleModelData>(
     compiler: SlangCompiler,
     hooks: impl FnOnce(&WebGpu, &SlangCompiler) -> Box<dyn MpmPipelineHooks<WebGpu, GpuModel>>,
     scene_builders: SceneBuilders<GpuModel>,
-    #[cfg(feature = "dim3")] up_axis: Vector3<f32>,
+    #[cfg(feature = "dim3")] up_axis: Vector,
 ) {
     run_with_hooks_and_ui(
         compiler,
@@ -373,7 +372,7 @@ pub async fn run_with_hooks_and_ui<GpuModel: GpuParticleModelData>(
         &SimulationStepResult,
         bool,
     ) -> Option<RunState>,
-    #[cfg(feature = "dim3")] up_axis: Vector3<f32>,
+    #[cfg(feature = "dim3")] up_axis: Vector,
 ) {
     let mut colliders_gfx = HashMap::new();
     let mut stage = Stage::new(compiler, hooks, scene_builders).await;
@@ -409,7 +408,9 @@ pub async fn run_with_hooks_and_ui<GpuModel: GpuParticleModelData>(
     };
     #[cfg(feature = "dim3")]
     {
-        camera3d.set_up_axis(up_axis);
+        camera3d.set_up_axis(kiss3d::nalgebra::Vector3::new(
+            up_axis.x, up_axis.y, up_axis.z,
+        ));
     }
     let mut camera2d = Sidescroll::new();
     #[cfg(feature = "dim2")]
@@ -632,17 +633,17 @@ fn update_colliders<GpuModel: GpuParticleModelData>(
             {
                 // TODO: here we are converting between nalgebra versions.
                 //       This can be simplified once kiss3d is updated to the latest nalgebra.
-                let tra = pose.translation.vector;
-                let rot = pose.rotation.into_inner();
+                let tra = pose.translation;
+                let rot = pose.rotation;
                 node.set_local_translation([tra.x, tra.y, tra.z].into());
                 node.set_local_rotation(kiss3d::nalgebra::Unit::new_unchecked(
-                    kiss3d::nalgebra::Quaternion::new(rot.w, rot.i, rot.j, rot.k),
+                    kiss3d::nalgebra::Quaternion::new(rot.w, rot.x, rot.y, rot.z),
                 ));
             }
             #[cfg(feature = "dim2")]
             {
-                let tra = pose.translation.vector;
-                let rot = pose.rotation.into_inner();
+                let tra = pose.translation;
+                let rot = pose.rotation;
                 node.set_local_translation([tra.x, tra.y].into());
                 node.set_local_rotation(kiss3d::nalgebra::Unit::new_unchecked(
                     kiss3d::nalgebra::Complex::new(rot.re, rot.im),
@@ -711,11 +712,11 @@ fn generate_collider_mesh(co_shape: &dyn Shape) -> Option<PlanarMesh> {
             for vox in voxels.voxels() {
                 if !vox.state.is_empty() {
                     let bid = vtx.len() as u32;
-                    let center = nalgebra::point![vox.center.x, vox.center.y];
-                    vtx.push(center + nalgebra::vector![sz.x, sz.y]);
-                    vtx.push(center + nalgebra::vector![-sz.x, sz.y]);
-                    vtx.push(center + nalgebra::vector![-sz.x, -sz.y]);
-                    vtx.push(center + nalgebra::vector![sz.x, -sz.y]);
+                    let center = vox.center;
+                    vtx.push(center + Vector::new(sz.x, sz.y));
+                    vtx.push(center + Vector::new(-sz.x, sz.y));
+                    vtx.push(center + Vector::new(-sz.x, -sz.y));
+                    vtx.push(center + Vector::new(sz.x, -sz.y));
                     idx.push([bid, bid + 1, bid + 2]);
                     idx.push([bid + 2, bid + 3, bid]);
                 }
@@ -753,7 +754,7 @@ fn generate_collider_mesh(co_shape: &dyn Shape) -> Option<PlanarMesh> {
 }
 
 #[cfg(feature = "dim2")]
-fn kiss3d_mesh_from_polyline(vertices: Vec<nalgebra::Point2<f32>>) -> PlanarMesh {
+fn kiss3d_mesh_from_polyline(vertices: Vec<Vector>) -> PlanarMesh {
     let n = vertices.len();
     let idx = (1..n as u32 - 1).map(|i| [0, i, i + 1]).collect();
     kiss3d_mesh((vertices, idx))
@@ -805,7 +806,7 @@ fn generate_collider_mesh(co_shape: &dyn Shape) -> Option<GpuMesh> {
 }
 
 #[cfg(feature = "dim3")]
-fn kiss3d_mesh(buffers: (Vec<nalgebra::Point3<f32>>, Vec<[u32; 3]>)) -> kiss3d::resource::GpuMesh {
+fn kiss3d_mesh(buffers: (Vec<Vector>, Vec<[u32; 3]>)) -> kiss3d::resource::GpuMesh {
     let (vtx, idx) = buffers;
     let kiss_vtx: Vec<_> = vtx
         .into_iter()
@@ -819,9 +820,7 @@ fn kiss3d_mesh(buffers: (Vec<nalgebra::Point3<f32>>, Vec<[u32; 3]>)) -> kiss3d::
 }
 
 #[cfg(feature = "dim2")]
-fn kiss3d_mesh(
-    buffers: (Vec<nalgebra::Point2<f32>>, Vec<[u32; 3]>),
-) -> kiss3d::resource::PlanarMesh {
+fn kiss3d_mesh(buffers: (Vec<Vector>, Vec<[u32; 3]>)) -> kiss3d::resource::PlanarMesh {
     let (vtx, idx) = buffers;
     let kiss_vtx: Vec<_> = vtx
         .into_iter()
diff --git a/src_testbed/prep_readback.rs b/src_testbed/prep_readback.rs
index 965c804..11a139d 100644
--- a/src_testbed/prep_readback.rs
+++ b/src_testbed/prep_readback.rs
@@ -1,5 +1,5 @@
 use bytemuck::{Pod, Zeroable};
-use nalgebra::Vector4;
+use glam::Vec4;
 use slang_hal::backend::{Backend, Encoder};
 use slang_hal::function::GpuFunction;
 use slang_hal::{BufferUsages, Shader, ShaderArgs};
@@ -15,7 +15,7 @@ use stensor::tensor::GpuTensor;
 #[derive(Default, Copy, Clone, Debug, Pod, Zeroable)]
 #[repr(C)]
 pub struct ReadbackData {
-    pub color: Vector4<f32>,
+    pub color: Vec4,
     pub deformation: nalgebra::Matrix2<f32>,
     pub position: nalgebra::Vector2<f32>,
     // NOTE: for now we are using explicit padding since
@@ -28,12 +28,12 @@ pub struct ReadbackData {
 #[derive(Default, Copy, Clone, Debug, Pod, Zeroable)]
 #[repr(C)]
 pub struct ReadbackData {
-    pub color: Vector4<f32>,
+    pub color: Vec4,
     // NOTE: for now we are using explicit padding since
     //       gpu buffer read based on Pod/bytemuck is much
     //       faster (about 20x) than with ShaderType/encase.
     pub deformation: nalgebra::Matrix4x3<f32>,
-    pub position: Vector4<f32>,
+    pub position: Vec4,
 }
 
 #[derive(Copy, Clone, PartialEq, Eq)]
@@ -104,7 +104,7 @@ pub struct PrepReadback<B: Backend> {
 
 pub struct GpuReadbackData<B: Backend> {
     pub mode: GpuTensor<RenderConfig, B>,
-    pub base_colors: GpuTensor<Vector4<f32>, B>,
+    pub base_colors: GpuTensor<Vec4, B>,
     pub instances: GpuTensor<ReadbackData, B>,
     pub instances_staging: GpuTensor<ReadbackData, B>,
 }
@@ -114,7 +114,7 @@ impl<B: Backend> GpuReadbackData<B> {
         backend: &B,
         num_particles: usize,
         mode: RenderMode,
-        custom_colors: Option<&[Vector4<f32>]>,
+        custom_colors: Option<&[Vec4]>,
     ) -> Result<Self, B::Error> {
         let config = mode.config();
         let palette = [
@@ -129,10 +129,10 @@ impl<B: Backend> GpuReadbackData<B> {
         let instances: Vec<_> = (0..num_particles)
             .map(|_| ReadbackData::default())
             .collect();
-        let base_colors: Vec<Vector4<f32>> = match custom_colors {
+        let base_colors: Vec<Vec4> = match custom_colors {
             Some(c) if c.len() == num_particles => c.to_vec(),
             _ => (0..num_particles)
-                .map(|i| palette[i % palette.len()].into())
+                .map(|i| Vec4::from_array(palette[i % palette.len()]))
                 .collect(),
         };
 
@@ -160,11 +160,11 @@ impl<B: Backend> GpuReadbackData<B> {
 #[derive(ShaderArgs)]
 struct PrepReadbackArgs<'a, B: Backend> {
     instances: &'a GpuTensor<ReadbackData, B>,
-    base_colors: &'a GpuTensor<Vector4<f32>, B>,
+    base_colors: &'a GpuTensor<Vec4, B>,
     particles_pos: &'a GpuTensor<ParticlePosition, B>,
     particles_kin: &'a GpuTensor<Kinematics, B>,
     particles_cdf: &'a GpuTensor<Cdf, B>,
-    particles_def_grad: &'a GpuTensor<Matrix<f32>, B>,
+    particles_def_grad: &'a GpuTensor<Matrix, B>,
     particles_props: &'a GpuTensor<ParticleProperties, B>,
     grid: &'a GpuTensor<GpuGridMetadata, B>,
     params: &'a GpuTensor<SimulationParams, B>,
diff --git a/src_testbed/step.rs b/src_testbed/step.rs
index 9f263f6..1e3c28b 100644
--- a/src_testbed/step.rs
+++ b/src_testbed/step.rs
@@ -4,14 +4,13 @@ use slang_hal::backend::{Backend, WebGpu};
 use slang_hal::BufferUsages;
 #[cfg(feature = "webgpu")]
 use slang_hal::GpuTimingResult;
-use slosh::rapier::na;
 use slosh::solver::{GpuParticleModelData, SimulationParams};
 use stensor::tensor::GpuTensor;
 
 /// Byte stride of one `Mat<f32>` element in the GPU def-grad buffer, matching
 /// the slang/WGSL structured-buffer layout. In 3D, `mat3x3<f32>` has three
 /// `vec4`-aligned columns for a total of 48 bytes (not the tightly-packed 36
-/// bytes of a CPU `nalgebra::Matrix3`). In 2D, `mat2x2<f32>` is already
+/// bytes of a CPU `glam::Mat3`). In 2D, `mat2x2<f32>` is already
 /// tightly packed at 16 bytes, so there is no mismatch.
 #[cfg(feature = "dim2")]
 pub const GPU_DEF_GRAD_STRIDE_BYTES: usize = 16;
@@ -287,9 +286,9 @@ impl<GpuModel: GpuParticleModelData> Stage<GpuModel> {
 
         // Copy deformation gradient buffer to staging for readback. We use the
         // GPU-side stride (see `GPU_DEF_GRAD_STRIDE_BYTES`) rather than
-        // `size_of::<Matrix<f32>>()`, which is only correct in 2D. In 3D the
+        // `size_of::<Matrix>()`, which is only correct in 2D. In 3D the
         // GPU mat3x3 has vec4-aligned columns (48 bytes), not the tightly
-        // packed 36 bytes of a CPU nalgebra Matrix3.
+        // packed 36 bytes of a CPU `glam::Mat3`.
         {
             let def_grad_buf = physics.data.particles.def_grad.buffer();
             let staging_buf = self.def_grad_staging.buffer();
@@ -337,7 +336,7 @@ impl<GpuModel: GpuParticleModelData> Stage<GpuModel> {
         // Step rapier to update kinematic bodies.
         let rapier = &mut self.physics.rapier_data;
         rapier.physics_pipeline.step(
-            &na::zero(),
+            slosh::math::Vector::ZERO,
             &rapier.params,
             &mut rapier.islands,
             &mut rapier.broad_phase,
diff --git a/validation_tests/Cargo.toml b/validation_tests/Cargo.toml
index 38839b0..355341f 100644
--- a/validation_tests/Cargo.toml
+++ b/validation_tests/Cargo.toml
@@ -18,6 +18,7 @@ webgpu = ["slosh3d/webgpu"]
 render = ["kiss3d"]
 
 [dependencies]
+glam = { workspace = true }
 nalgebra = { workspace = true }
 wgpu = { workspace = true }
 slang-hal = { workspace = true }
@@ -38,7 +39,10 @@ plotters = "0.3"
 rstar = "0.12"
 
 # Rendering dependencies (optional)
-kiss3d = { version = "0.37", features = ["vertex_index_u32", "egui"], optional = true }
+kiss3d = { version = "0.37", features = [
+    "vertex_index_u32",
+    "egui",
+], optional = true }
 
 [dependencies.pollster]
 version = "0.4"