Skip to main content

oxiphysics_gpu/deformable_gpu/
functions.rs

1//! Auto-generated module
2//!
3//! 🤖 Generated with [SplitRS](https://github.com/cool-japan/splitrs)
4
5#[cfg(test)]
6use super::types::*;
7
8#[inline]
9pub(super) fn dot3(a: [f64; 3], b: [f64; 3]) -> f64 {
10    a[0] * b[0] + a[1] * b[1] + a[2] * b[2]
11}
12#[inline]
13pub(super) fn cross3(a: [f64; 3], b: [f64; 3]) -> [f64; 3] {
14    [
15        a[1] * b[2] - a[2] * b[1],
16        a[2] * b[0] - a[0] * b[2],
17        a[0] * b[1] - a[1] * b[0],
18    ]
19}
20#[inline]
21pub(super) fn length3(v: [f64; 3]) -> f64 {
22    dot3(v, v).sqrt()
23}
24#[inline]
25pub(super) fn normalize3(v: [f64; 3]) -> [f64; 3] {
26    let len = length3(v);
27    if len < 1e-30 {
28        return [0.0; 3];
29    }
30    scale3(v, 1.0 / len)
31}
32#[inline]
33pub(super) fn sub3(a: [f64; 3], b: [f64; 3]) -> [f64; 3] {
34    [a[0] - b[0], a[1] - b[1], a[2] - b[2]]
35}
36#[inline]
37pub(super) fn add3(a: [f64; 3], b: [f64; 3]) -> [f64; 3] {
38    [a[0] + b[0], a[1] + b[1], a[2] + b[2]]
39}
40#[inline]
41pub(super) fn scale3(v: [f64; 3], s: f64) -> [f64; 3] {
42    [v[0] * s, v[1] * s, v[2] * s]
43}
44/// Compute the squared length of a 3-D vector.
45#[inline]
46pub(super) fn len2_3(v: [f64; 3]) -> f64 {
47    dot3(v, v)
48}
49/// Rotate a point `p` by unit quaternion `q = [w, x, y, z]`.
50pub(super) fn quat_rotate(q: [f64; 4], p: [f64; 3]) -> [f64; 3] {
51    let [w, x, y, z] = q;
52    let qv = [x, y, z];
53    let t = scale3(cross3(qv, p), 2.0);
54    add3(add3(p, scale3(t, w)), cross3(qv, t))
55}
56#[cfg(test)]
57mod tests {
58    use super::*;
59    fn make_simple_mesh(n: usize) -> DeformableGpuMesh {
60        let mut mesh = DeformableGpuMesh::new();
61        for i in 0..n {
62            let x = i as f64;
63            mesh.add_vertex(DeformableVertex::new([x, 0.0, 0.0]));
64        }
65        mesh
66    }
67    #[test]
68    fn test_mesh_add_vertex() {
69        let mut mesh = DeformableGpuMesh::new();
70        let idx = mesh.add_vertex(DeformableVertex::new([1.0, 2.0, 3.0]));
71        assert_eq!(idx, 0);
72        assert_eq!(mesh.vertices.len(), 1);
73    }
74    #[test]
75    fn test_mesh_add_face() {
76        let mut mesh = make_simple_mesh(4);
77        mesh.add_face(0, 1, 2);
78        assert_eq!(mesh.faces.len(), 1);
79        assert_eq!(mesh.faces[0].a, 0);
80    }
81    #[test]
82    fn test_mesh_aabb() {
83        let mut mesh = DeformableGpuMesh::new();
84        mesh.add_vertex(DeformableVertex::new([0.0, 0.0, 0.0]));
85        mesh.add_vertex(DeformableVertex::new([1.0, 2.0, 3.0]));
86        let (lo, hi) = mesh.aabb();
87        assert!((lo[0]).abs() < 1e-10);
88        assert!((hi[2] - 3.0).abs() < 1e-10);
89    }
90    #[test]
91    fn test_mesh_aabb_single_vertex() {
92        let mut mesh = DeformableGpuMesh::new();
93        mesh.add_vertex(DeformableVertex::new([5.0, 6.0, 7.0]));
94        let (lo, hi) = mesh.aabb();
95        assert!((lo[0] - 5.0).abs() < 1e-10);
96        assert!((hi[0] - 5.0).abs() < 1e-10);
97    }
98    #[test]
99    fn test_mesh_face_normals() {
100        let mut mesh = DeformableGpuMesh::new();
101        mesh.add_vertex(DeformableVertex::new([0.0, 0.0, 0.0]));
102        mesh.add_vertex(DeformableVertex::new([1.0, 0.0, 0.0]));
103        mesh.add_vertex(DeformableVertex::new([0.0, 1.0, 0.0]));
104        mesh.add_face(0, 1, 2);
105        let normals = mesh.compute_face_normals();
106        assert_eq!(normals.len(), 1);
107        assert!((normals[0][2] - 1.0).abs() < 1e-10);
108    }
109    #[test]
110    fn test_blend_shape_apply() {
111        let mut mesh = DeformableGpuMesh::new();
112        mesh.add_vertex(DeformableVertex::new([0.0, 0.0, 0.0]));
113        let bs = BlendShape::new("test", vec![[0.0, 1.0, 0.0]]);
114        let idx = mesh.add_blend_shape(bs);
115        mesh.blend_shapes[idx].weight = 1.0;
116        mesh.apply_blend_shapes();
117        assert!((mesh.vertices[0].position[1] - 1.0).abs() < 1e-10);
118    }
119    #[test]
120    fn test_blend_shape_zero_weight() {
121        let mut mesh = DeformableGpuMesh::new();
122        mesh.add_vertex(DeformableVertex::new([0.0, 0.0, 0.0]));
123        let bs = BlendShape::new("test", vec![[0.0, 1.0, 0.0]]);
124        mesh.add_blend_shape(bs);
125        mesh.apply_blend_shapes();
126        assert!((mesh.vertices[0].position[1]).abs() < 1e-10);
127    }
128    #[test]
129    fn test_joint_transform_identity() {
130        let jt = JointTransform::identity();
131        let p = [1.0, 2.0, 3.0];
132        let tp = jt.transform_point(p);
133        for k in 0..3 {
134            assert!((tp[k] - p[k]).abs() < 1e-10, "k={k}");
135        }
136    }
137    #[test]
138    fn test_joint_transform_translation() {
139        let mut jt = JointTransform::identity();
140        jt.row0[3] = 5.0;
141        let tp = jt.transform_point([0.0, 0.0, 0.0]);
142        assert!((tp[0] - 5.0).abs() < 1e-10);
143    }
144    #[test]
145    fn test_dualquat_identity_transform() {
146        let dq = DualQuat::identity();
147        let p = [1.0, 2.0, 3.0];
148        let tp = dq.transform_point(p);
149        for k in 0..3 {
150            assert!((tp[k] - p[k]).abs() < 1e-9, "k={k}, got {}", tp[k]);
151        }
152    }
153    #[test]
154    fn test_dualquat_pure_translation() {
155        let dq = DualQuat::from_translation([3.0, 0.0, 0.0]);
156        let p = [0.0, 0.0, 0.0];
157        let tp = dq.transform_point(p);
158        assert!((tp[0] - 3.0).abs() < 1e-9, "got {}", tp[0]);
159    }
160    #[test]
161    fn test_dualquat_normalize_identity() {
162        let dq = DualQuat::identity();
163        let dn = dq.normalize();
164        assert!((dn.real[0] - 1.0).abs() < 1e-10);
165    }
166    #[test]
167    fn test_skinning_lbs_identity() {
168        let kernel = SkinningKernel::new(SkinningMode::LinearBlend, 1);
169        let v = DeformableVertex::new([1.0, 2.0, 3.0]);
170        let result = kernel.apply_lbs(&v);
171        for (result_val, rest_val) in result.iter().zip(v.rest_position.iter()) {
172            assert!((result_val - rest_val).abs() < 1e-10);
173        }
174    }
175    #[test]
176    fn test_skinning_dqs_identity() {
177        let kernel = SkinningKernel::new(SkinningMode::DualQuaternion, 1);
178        let v = DeformableVertex::new([1.0, 2.0, 3.0]);
179        let result = kernel.apply_dqs(&v);
180        for (k, (result_val, rest_val)) in result.iter().zip(v.rest_position.iter()).enumerate() {
181            assert!((result_val - rest_val).abs() < 1e-9, "k={k}");
182        }
183    }
184    #[test]
185    fn test_skinning_dispatch_lbs() {
186        let mut kernel = SkinningKernel::new(SkinningMode::LinearBlend, 1);
187        kernel.joints[0].row1[3] = 5.0;
188        let mut mesh = make_simple_mesh(3);
189        kernel.dispatch(&mut mesh);
190        for v in &mesh.vertices {
191            assert!((v.position[1] - 5.0).abs() < 1e-10);
192        }
193    }
194    #[test]
195    fn test_skinning_dispatch_dqs() {
196        let mut kernel = SkinningKernel::new(SkinningMode::DualQuaternion, 1);
197        kernel.dq_joints[0] = DualQuat::from_translation([0.0, 5.0, 0.0]);
198        let mut mesh = make_simple_mesh(3);
199        kernel.dispatch(&mut mesh);
200        for v in &mesh.vertices {
201            assert!((v.position[1] - 5.0).abs() < 1e-9, "got {}", v.position[1]);
202        }
203    }
204    #[test]
205    fn test_xpbd_no_constraint() {
206        let kernel = XPBDGpuKernel::new(4);
207        let mut pos = vec![[0.0f64; 3], [1.0, 0.0, 0.0]];
208        let inv_m = vec![1.0, 1.0];
209        kernel.solve_step(&mut pos, &inv_m, 0.01);
210        assert!((pos[0][0]).abs() < 1e-10);
211        assert!((pos[1][0] - 1.0).abs() < 1e-10);
212    }
213    #[test]
214    fn test_xpbd_distance_constraint_satisfied() {
215        let mut kernel = XPBDGpuKernel::new(20);
216        kernel.add_constraint(XpbdDistConstraint::new(0, 1, 1.0, 0.0));
217        let mut pos = vec![[0.0f64; 3], [2.0, 0.0, 0.0]];
218        let inv_m = vec![1.0, 1.0];
219        kernel.solve_step(&mut pos, &inv_m, 0.01);
220        let dist = length3(sub3(pos[1], pos[0]));
221        assert!((dist - 1.0).abs() < 1e-6, "dist = {dist}");
222    }
223    #[test]
224    fn test_xpbd_fixed_particle() {
225        let mut kernel = XPBDGpuKernel::new(10);
226        kernel.add_constraint(XpbdDistConstraint::new(0, 1, 1.0, 0.0));
227        let mut pos = vec![[0.0f64; 3], [5.0, 0.0, 0.0]];
228        let inv_m = vec![0.0, 1.0];
229        kernel.solve_step(&mut pos, &inv_m, 0.01);
230        assert!((pos[0][0]).abs() < 1e-10);
231        let dist = length3(sub3(pos[1], pos[0]));
232        assert!((dist - 1.0).abs() < 1e-6, "dist = {dist}");
233    }
234    #[test]
235    fn test_xpbd_integrate_gravity() {
236        let kernel = XPBDGpuKernel::new(1);
237        let mut pos = vec![[0.0f64; 3]];
238        let mut vel = vec![[0.0f64; 3]];
239        let inv_m = vec![1.0];
240        let gravity = [0.0, -9.81, 0.0];
241        let dt = 0.1;
242        kernel.integrate(&mut pos, &mut vel, &inv_m, gravity, dt);
243        assert!((vel[0][1] - (-9.81 * dt)).abs() < 1e-10);
244        assert!((pos[0][1] - (-9.81 * dt * dt)).abs() < 1e-10);
245    }
246    #[test]
247    fn test_xpbd_integrate_fixed_particle() {
248        let kernel = XPBDGpuKernel::new(1);
249        let mut pos = vec![[0.0f64; 3]];
250        let mut vel = vec![[0.0f64; 3]];
251        let inv_m = vec![0.0];
252        kernel.integrate(&mut pos, &mut vel, &inv_m, [0.0, -9.81, 0.0], 0.1);
253        assert!((pos[0][1]).abs() < 1e-10);
254    }
255    #[test]
256    fn test_fem_forces_at_rest() {
257        let mut fem = FEMGpuKernel::new(4);
258        let elem = TetElement::new([0, 1, 2, 3], 1.0 / 6.0, 1e6, 0.3);
259        fem.add_element(elem);
260        let rest = vec![
261            [0.0f64, 0.0, 0.0],
262            [1.0, 0.0, 0.0],
263            [0.0, 1.0, 0.0],
264            [0.0, 0.0, 1.0],
265        ];
266        let forces = fem.compute_forces(&rest, &rest);
267        for f in &forces {
268            let mag = length3(*f);
269            assert!(mag < 1e-10, "mag = {mag}");
270        }
271    }
272    #[test]
273    fn test_fem_forces_nonzero_when_stretched() {
274        let mut fem = FEMGpuKernel::new(4);
275        let elem = TetElement::new([0, 1, 2, 3], 1.0 / 6.0, 1e6, 0.3);
276        fem.add_element(elem);
277        let rest = vec![
278            [0.0f64, 0.0, 0.0],
279            [1.0, 0.0, 0.0],
280            [0.0, 1.0, 0.0],
281            [0.0, 0.0, 1.0],
282        ];
283        let mut deformed = rest.clone();
284        deformed[1] = [2.0, 0.0, 0.0];
285        let forces = fem.compute_forces(&deformed, &rest);
286        let total_mag: f64 = forces.iter().map(|f| length3(*f)).sum();
287        assert!(total_mag > 1e-6, "expected nonzero forces");
288    }
289    #[test]
290    fn test_fem_strain_energy_zero_at_rest() {
291        let mut fem = FEMGpuKernel::new(4);
292        let elem = TetElement::new([0, 1, 2, 3], 1.0, 1e5, 0.3);
293        fem.add_element(elem);
294        let rest = vec![
295            [0.0f64, 0.0, 0.0],
296            [1.0, 0.0, 0.0],
297            [0.0, 1.0, 0.0],
298            [0.0, 0.0, 1.0],
299        ];
300        let energy = fem.compute_strain_energy(&rest, &rest);
301        assert!(energy < 1e-10, "energy = {energy}");
302    }
303    #[test]
304    fn test_fem_lame_params() {
305        let elem = TetElement::new([0, 1, 2, 3], 1.0, 1e6, 0.25);
306        let (lambda, mu) = elem.lame_params();
307        assert!(lambda > 0.0);
308        assert!(mu > 0.0);
309    }
310    #[test]
311    fn test_collision_signed_distance_above() {
312        let cr = CollisionResponseGpu::new(CollisionResponseParams::default());
313        let sd = cr.signed_distance_to_triangle(
314            [0.0, 1.0, 0.0],
315            [0.0, 0.0, 0.0],
316            [1.0, 0.0, 0.0],
317            [0.0, 0.0, 1.0],
318        );
319        assert!(sd > 0.0, "sd = {sd}");
320    }
321    #[test]
322    fn test_collision_signed_distance_below() {
323        let cr = CollisionResponseGpu::new(CollisionResponseParams::default());
324        let sd = cr.signed_distance_to_triangle(
325            [0.0, -0.5, 0.0],
326            [0.0, 0.0, 0.0],
327            [1.0, 0.0, 0.0],
328            [0.0, 0.0, 1.0],
329        );
330        assert!(sd < 0.0, "sd = {sd}");
331    }
332    #[test]
333    fn test_collision_penalty_forces_below_threshold() {
334        let cr = CollisionResponseGpu::new(CollisionResponseParams::new(1e4, 10.0, 1.0));
335        let positions = vec![[0.5, -0.001, 0.5f64]];
336        let velocities = vec![[0.0f64; 3]];
337        let faces = vec![([0.0, 0.0, 0.0f64], [1.0, 0.0, 0.0f64], [0.0, 0.0, 1.0f64])];
338        let mut forces = vec![[0.0f64; 3]];
339        cr.apply_penalty_forces(&positions, &velocities, &faces, &mut forces);
340        assert!(forces[0][1] > 0.0, "expected upward penalty force");
341    }
342    #[test]
343    fn test_collision_penalty_forces_above_plane() {
344        let cr = CollisionResponseGpu::new(CollisionResponseParams::default());
345        let positions = vec![[0.5, 1.0, 0.5f64]];
346        let velocities = vec![[0.0f64; 3]];
347        let faces = vec![([0.0, 0.0, 0.0f64], [1.0, 0.0, 0.0f64], [0.0, 0.0, 1.0f64])];
348        let mut forces = vec![[0.0f64; 3]];
349        cr.apply_penalty_forces(&positions, &velocities, &faces, &mut forces);
350        let mag = length3(forces[0]);
351        assert!(mag < 1e-10, "expected no force, got mag = {mag}");
352    }
353    #[test]
354    fn test_collision_sphere_penalty_inside() {
355        let cr = CollisionResponseGpu::new(CollisionResponseParams::new(1e4, 0.0, 0.01));
356        let positions = vec![[0.0, 0.5, 0.0f64]];
357        let velocities = vec![[0.0f64; 3]];
358        let mut forces = vec![[0.0f64; 3]];
359        cr.apply_sphere_penalty(&positions, &velocities, [0.0, 0.0, 0.0], 1.0, &mut forces);
360        let mag = length3(forces[0]);
361        assert!(mag > 0.0, "expected sphere penalty force");
362    }
363    #[test]
364    fn test_collision_sphere_penalty_outside() {
365        let cr = CollisionResponseGpu::new(CollisionResponseParams::default());
366        let positions = vec![[0.0, 2.0, 0.0f64]];
367        let velocities = vec![[0.0f64; 3]];
368        let mut forces = vec![[0.0f64; 3]];
369        cr.apply_sphere_penalty(&positions, &velocities, [0.0, 0.0, 0.0], 1.0, &mut forces);
370        let mag = length3(forces[0]);
371        assert!(mag < 1e-10, "expected no sphere force");
372    }
373    #[test]
374    fn test_pipeline_creation() {
375        let mesh = make_simple_mesh(4);
376        let config = PhysicsPipelineConfig::default();
377        let pipeline = PhysicsGpuPipeline::new(config, mesh);
378        assert_eq!(pipeline.mesh.vertices.len(), 4);
379        assert!((pipeline.time).abs() < 1e-10);
380    }
381    #[test]
382    fn test_pipeline_step_advances_time() {
383        let mesh = make_simple_mesh(4);
384        let config = PhysicsPipelineConfig::default();
385        let mut pipeline = PhysicsGpuPipeline::new(config, mesh);
386        pipeline.step();
387        assert!(pipeline.time > 0.0);
388    }
389    #[test]
390    fn test_pipeline_gravity_falls() {
391        let mut mesh = DeformableGpuMesh::new();
392        let v = DeformableVertex::new([0.0, 10.0, 0.0]);
393        mesh.add_vertex(v);
394        let config = PhysicsPipelineConfig {
395            dt: 0.1,
396            substeps: 1,
397            xpbd_enabled: true,
398            fem_enabled: false,
399            collision_enabled: false,
400            ..Default::default()
401        };
402        let mut pipeline = PhysicsGpuPipeline::new(config, mesh);
403        pipeline.step();
404        assert!(
405            pipeline.mesh.vertices[0].position[1] < 10.0,
406            "y = {}",
407            pipeline.mesh.vertices[0].position[1]
408        );
409    }
410    #[test]
411    fn test_pipeline_centre_of_mass() {
412        let mut mesh = DeformableGpuMesh::new();
413        mesh.add_vertex(DeformableVertex::new([0.0, 0.0, 0.0]));
414        mesh.add_vertex(DeformableVertex::new([2.0, 0.0, 0.0]));
415        let config = PhysicsPipelineConfig::default();
416        let pipeline = PhysicsGpuPipeline::new(config, mesh);
417        let com = pipeline.centre_of_mass();
418        assert!((com[0] - 1.0).abs() < 1e-10);
419    }
420    #[test]
421    fn test_pipeline_kinetic_energy_zero_at_rest() {
422        let mesh = make_simple_mesh(4);
423        let config = PhysicsPipelineConfig::default();
424        let pipeline = PhysicsGpuPipeline::new(config, mesh);
425        assert!((pipeline.kinetic_energy()).abs() < 1e-10);
426    }
427    #[test]
428    fn test_pipeline_xpbd_constraint_spring() {
429        let mut mesh = DeformableGpuMesh::new();
430        let mut v0 = DeformableVertex::new([0.0, 0.0, 0.0]);
431        v0.mass = 1e30;
432        mesh.add_vertex(v0);
433        let v1 = DeformableVertex::new([5.0, 0.0, 0.0]);
434        mesh.add_vertex(v1);
435        let config = PhysicsPipelineConfig {
436            gravity: [0.0; 3],
437            dt: 0.01,
438            substeps: 10,
439            collision_enabled: false,
440            fem_enabled: false,
441            xpbd_enabled: true,
442        };
443        let mut pipeline = PhysicsGpuPipeline::new(config, mesh);
444        pipeline
445            .xpbd
446            .add_constraint(XpbdDistConstraint::new(0, 1, 1.0, 0.0));
447        pipeline.step();
448        let dx = pipeline.mesh.vertices[1].position[0];
449        assert!(dx < 5.0, "particle should move, dx = {dx}");
450    }
451    #[test]
452    fn test_xpbd_multiple_constraints() {
453        let mut kernel = XPBDGpuKernel::new(50);
454        kernel.add_constraint(XpbdDistConstraint::new(0, 1, 1.0, 0.0));
455        kernel.add_constraint(XpbdDistConstraint::new(1, 2, 1.0, 0.0));
456        let mut pos = vec![[0.0f64; 3], [3.0, 0.0, 0.0], [6.0, 0.0, 0.0]];
457        let inv_m = vec![0.0, 1.0, 1.0];
458        kernel.solve_step(&mut pos, &inv_m, 0.01);
459        let d01 = length3(sub3(pos[1], pos[0]));
460        let d12 = length3(sub3(pos[2], pos[1]));
461        assert!((d01 - 1.0).abs() < 1e-5, "d01 = {d01}");
462        assert!((d12 - 1.0).abs() < 1e-5, "d12 = {d12}");
463    }
464    #[test]
465    fn test_blend_shape_partial_weight() {
466        let mut mesh = DeformableGpuMesh::new();
467        mesh.add_vertex(DeformableVertex::new([0.0, 0.0, 0.0]));
468        let bs = BlendShape::new("half", vec![[0.0, 2.0, 0.0]]);
469        let idx = mesh.add_blend_shape(bs);
470        mesh.blend_shapes[idx].weight = 0.5;
471        mesh.apply_blend_shapes();
472        assert!((mesh.vertices[0].position[1] - 1.0).abs() < 1e-10);
473    }
474    #[test]
475    fn test_fem_add_element_updates_num_nodes() {
476        let mut fem = FEMGpuKernel::new(0);
477        fem.add_element(TetElement::new([0, 1, 2, 10], 1.0, 1e5, 0.3));
478        assert!(fem.num_nodes >= 11);
479    }
480    #[test]
481    fn test_collision_params_default() {
482        let p = CollisionResponseParams::default();
483        assert!(p.penalty_stiffness > 0.0);
484        assert!(p.damping >= 0.0);
485        assert!(p.contact_threshold > 0.0);
486    }
487    #[test]
488    fn test_deformable_vertex_default_weights() {
489        let v = DeformableVertex::new([1.0, 2.0, 3.0]);
490        assert!((v.joint_weights[0] - 1.0).abs() < 1e-10);
491        assert!((v.joint_weights[1]).abs() < 1e-10);
492    }
493    #[test]
494    fn test_quat_rotate_identity() {
495        let q = [1.0, 0.0, 0.0, 0.0];
496        let p = [3.0, 4.0, 5.0];
497        let r = quat_rotate(q, p);
498        for k in 0..3 {
499            assert!((r[k] - p[k]).abs() < 1e-9);
500        }
501    }
502    #[test]
503    fn test_skinning_mode_eq() {
504        assert_eq!(SkinningMode::LinearBlend, SkinningMode::LinearBlend);
505        assert_ne!(SkinningMode::LinearBlend, SkinningMode::DualQuaternion);
506    }
507    #[test]
508    fn test_pipeline_multiple_steps() {
509        let mesh = make_simple_mesh(2);
510        let config = PhysicsPipelineConfig {
511            gravity: [0.0, -9.81, 0.0],
512            dt: 0.016,
513            substeps: 2,
514            collision_enabled: false,
515            fem_enabled: false,
516            xpbd_enabled: true,
517        };
518        let mut pipeline = PhysicsGpuPipeline::new(config, mesh);
519        for _ in 0..10 {
520            pipeline.step();
521        }
522        assert!((pipeline.time - 0.16).abs() < 1e-6);
523    }
524    #[test]
525    fn test_gpu_node_new() {
526        let node = DeformableGpuNode::new([1.0, 2.0, 3.0], 2.5_f64);
527        assert!((node.mass - 2.5_f64).abs() < 1e-10);
528        assert_eq!(node.vel, [0.0; 3]);
529        assert_eq!(node.force, [0.0; 3]);
530    }
531    #[test]
532    fn test_gpu_node_apply_force() {
533        let mut node = DeformableGpuNode::new([0.0; 3], 1.0_f64);
534        node.apply_force([1.0, 0.0, 0.0]);
535        node.apply_force([0.0, 2.0, 0.0]);
536        assert!((node.force[0] - 1.0_f64).abs() < 1e-10);
537        assert!((node.force[1] - 2.0_f64).abs() < 1e-10);
538    }
539    #[test]
540    fn test_gpu_node_reset_force() {
541        let mut node = DeformableGpuNode::new([0.0; 3], 1.0_f64);
542        node.apply_force([5.0, 5.0, 5.0]);
543        node.reset_force();
544        assert_eq!(node.force, [0.0; 3]);
545    }
546    #[test]
547    fn test_gpu_node_kinetic_energy() {
548        let mut node = DeformableGpuNode::new([0.0; 3], 2.0_f64);
549        node.vel = [3.0, 4.0, 0.0];
550        assert!((node.kinetic_energy() - 25.0_f64).abs() < 1e-10);
551    }
552    #[test]
553    fn test_gpu_buffer_push_len() {
554        let mut buf = DeformableGpuBuffer::new();
555        buf.push([0.0; 3], [0.0; 3], 1.0_f64);
556        buf.push([1.0, 0.0, 0.0], [0.0; 3], 1.0_f64);
557        assert_eq!(buf.len(), 2);
558        assert!(!buf.is_empty());
559    }
560    #[test]
561    fn test_gpu_buffer_from_nodes() {
562        let nodes = vec![
563            DeformableGpuNode::new([0.0; 3], 1.0_f64),
564            DeformableGpuNode::new([1.0, 0.0, 0.0], 2.0_f64),
565        ];
566        let buf = DeformableGpuBuffer::from_nodes(&nodes);
567        assert_eq!(buf.len(), 2);
568        assert!((buf.masses[1] - 2.0_f64).abs() < 1e-10);
569        assert!((buf.inv_masses[1] - 0.5_f64).abs() < 1e-10);
570    }
571    #[test]
572    fn test_gpu_buffer_reset_forces() {
573        let mut buf = DeformableGpuBuffer::new();
574        buf.push([0.0; 3], [0.0; 3], 1.0_f64);
575        buf.forces[0] = [5.0, 5.0, 5.0];
576        buf.reset_forces();
577        assert_eq!(buf.forces[0], [0.0; 3]);
578    }
579    #[test]
580    fn test_gpu_buffer_apply_gravity() {
581        let mut buf = DeformableGpuBuffer::new();
582        buf.push([0.0; 3], [0.0; 3], 2.0_f64);
583        buf.apply_gravity([0.0, -9.81_f64, 0.0]);
584        assert!((buf.forces[0][1] - (-19.62_f64)).abs() < 1e-5);
585    }
586    #[test]
587    fn test_gpu_buffer_integrate() {
588        let mut buf = DeformableGpuBuffer::new();
589        buf.push([0.0; 3], [0.0; 3], 1.0_f64);
590        buf.forces[0] = [0.0, -9.81_f64, 0.0];
591        buf.integrate(0.1_f64);
592        assert!((buf.velocities[0][1] - (-0.981_f64)).abs() < 1e-6);
593    }
594    #[test]
595    fn test_gpu_buffer_total_ke() {
596        let mut buf = DeformableGpuBuffer::new();
597        buf.push([0.0; 3], [1.0, 0.0, 0.0], 2.0_f64);
598        assert!((buf.total_kinetic_energy() - 1.0_f64).abs() < 1e-10);
599    }
600    #[test]
601    fn test_gpu_buffer_flatten_positions() {
602        let mut buf = DeformableGpuBuffer::new();
603        buf.push([1.0, 2.0, 3.0], [0.0; 3], 1.0_f64);
604        let flat = buf.flatten_positions();
605        assert_eq!(flat, vec![1.0, 2.0, 3.0]);
606    }
607    #[test]
608    fn test_gpu_buffer_fixed_particle_inv_mass() {
609        let mut buf = DeformableGpuBuffer::new();
610        buf.push([0.0; 3], [0.0; 3], 0.0_f64);
611        assert!((buf.inv_masses[0]).abs() < 1e-15);
612    }
613    #[test]
614    fn test_xpbd_solver_distance_constraint() {
615        let mut solver = XPBDGpuSolver::new(20);
616        solver.add_distance(0, 1, 1.0_f64, 0.0_f64);
617        let mut buf = DeformableGpuBuffer::new();
618        buf.push([0.0; 3], [0.0; 3], 1.0_f64);
619        buf.push([3.0, 0.0, 0.0], [0.0; 3], 1.0_f64);
620        solver.solve(&mut buf, 0.01_f64);
621        let dist = length3(sub3(buf.positions[1], buf.positions[0]));
622        assert!((dist - 1.0_f64).abs() < 0.01_f64, "dist = {dist}");
623    }
624    #[test]
625    fn test_xpbd_solver_no_constraint_unchanged() {
626        let solver = XPBDGpuSolver::new(5);
627        let mut buf = DeformableGpuBuffer::new();
628        buf.push([0.0; 3], [0.0; 3], 1.0_f64);
629        buf.push([2.0, 0.0, 0.0], [0.0; 3], 1.0_f64);
630        solver.solve(&mut buf, 0.01_f64);
631        assert!((buf.positions[0][0]).abs() < 1e-10);
632        assert!((buf.positions[1][0] - 2.0_f64).abs() < 1e-10);
633    }
634    #[test]
635    fn test_xpbd_solver_shape_matching() {
636        let mut solver = XPBDGpuSolver::new(20);
637        let rest = vec![[0.0; 3], [1.0, 0.0, 0.0]];
638        solver.add_shape(ShapeMatchConstraint::new(vec![0, 1], rest, 1.0_f64));
639        let mut buf = DeformableGpuBuffer::new();
640        buf.push([0.0; 3], [0.0; 3], 1.0_f64);
641        buf.push([3.0, 0.0, 0.0], [0.0; 3], 1.0_f64);
642        solver.solve(&mut buf, 0.01_f64);
643        let gap = (buf.positions[1][0] - buf.positions[0][0]).abs();
644        assert!(gap < 3.0_f64, "gap = {gap}");
645    }
646    #[test]
647    fn test_fem_solver_rest_zero_energy() {
648        let rest = [
649            [0.0, 0.0, 0.0f64],
650            [1.0, 0.0, 0.0],
651            [0.0, 1.0, 0.0],
652            [0.0, 0.0, 1.0],
653        ];
654        let elem = CorotationalElement::from_rest([0, 1, 2, 3], &rest, 1e6_f64, 0.3_f64);
655        let mut solver = FEMGpuSolver::new(4);
656        solver.add_element(elem);
657        let mut buf = DeformableGpuBuffer::new();
658        for &p in &rest {
659            buf.push(p, [0.0; 3], 1.0_f64);
660        }
661        let energy = solver.strain_energy(&buf);
662        assert!(energy < 1e-5_f64, "energy = {energy}");
663    }
664    #[test]
665    fn test_fem_solver_assemble_forces_at_rest() {
666        let rest = [
667            [0.0, 0.0, 0.0f64],
668            [1.0, 0.0, 0.0],
669            [0.0, 1.0, 0.0],
670            [0.0, 0.0, 1.0],
671        ];
672        let elem = CorotationalElement::from_rest([0, 1, 2, 3], &rest, 1e5_f64, 0.3_f64);
673        let mut solver = FEMGpuSolver::new(4);
674        solver.add_element(elem);
675        let mut buf = DeformableGpuBuffer::new();
676        for &p in &rest {
677            buf.push(p, [0.0; 3], 1.0_f64);
678        }
679        buf.reset_forces();
680        solver.assemble_forces(&mut buf);
681        for i in 0..4 {
682            let mag = length3(buf.forces[i]);
683            assert!(mag < 1e-5_f64, "node {i} force mag = {mag}");
684        }
685    }
686    #[test]
687    fn test_fem_corotational_element_lame() {
688        let rest = [[0.0; 3], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]];
689        let elem = CorotationalElement::from_rest([0, 1, 2, 3], &rest, 1e6_f64, 0.25_f64);
690        let (lambda, mu) = elem.lame();
691        assert!(lambda > 0.0_f64);
692        assert!(mu > 0.0_f64);
693    }
694    #[test]
695    fn test_fem_solver_num_nodes_update() {
696        let rest = [[0.0; 3]; 4];
697        let elem = CorotationalElement::from_rest([0, 1, 9, 3], &rest, 1e5_f64, 0.3_f64);
698        let mut solver = FEMGpuSolver::new(0);
699        solver.add_element(elem);
700        assert!(solver.num_nodes >= 10);
701    }
702    #[test]
703    fn test_gpu_collision_plane_detect() {
704        let cr = GPUCollisionResponse::default_params();
705        let positions = vec![[0.0, -0.01_f64, 0.0]];
706        let contacts = cr.detect_plane_contacts(&positions, [0.0, 1.0, 0.0], 0.0_f64);
707        assert!(!contacts.is_empty());
708        assert_eq!(contacts[0].particle_idx, 0);
709    }
710    #[test]
711    fn test_gpu_collision_plane_no_contact_above() {
712        let cr = GPUCollisionResponse::default_params();
713        let positions = vec![[0.0, 1.0_f64, 0.0]];
714        let contacts = cr.detect_plane_contacts(&positions, [0.0, 1.0, 0.0], 0.0_f64);
715        assert!(contacts.is_empty());
716    }
717    #[test]
718    fn test_gpu_collision_sphere_detect() {
719        let cr = GPUCollisionResponse::default_params();
720        let positions = vec![[0.5_f64, 0.0, 0.0]];
721        let contacts = cr.detect_sphere_contacts(&positions, [0.0; 3], 1.0_f64);
722        assert!(!contacts.is_empty());
723    }
724    #[test]
725    fn test_gpu_collision_resolve_positions() {
726        let cr = GPUCollisionResponse::default_params();
727        let mut buf = DeformableGpuBuffer::new();
728        buf.push([0.0, -0.05_f64, 0.0], [0.0; 3], 1.0_f64);
729        let contacts = cr.detect_plane_contacts(&buf.positions, [0.0, 1.0, 0.0], 0.0_f64);
730        cr.resolve_positions(&mut buf, &contacts);
731        assert!(
732            buf.positions[0][1] >= 0.0_f64 - 1e-10,
733            "y = {}",
734            buf.positions[0][1]
735        );
736    }
737    #[test]
738    fn test_gpu_collision_resolve_velocities() {
739        let cr = GPUCollisionResponse::new(1e4_f64, 10.0_f64, 0.01_f64, 0.5_f64);
740        let mut buf = DeformableGpuBuffer::new();
741        buf.push([0.0, -0.001_f64, 0.0], [0.0, -2.0_f64, 0.0], 1.0_f64);
742        let contacts = cr.detect_plane_contacts(&buf.positions, [0.0, 1.0, 0.0], 0.0_f64);
743        cr.resolve_velocities(&mut buf, &contacts);
744        assert!(
745            buf.velocities[0][1] > 0.0_f64,
746            "vel = {}",
747            buf.velocities[0][1]
748        );
749    }
750    #[test]
751    fn test_deformable_pipeline_creation() {
752        let pipeline = DeformableGpuPipeline::new([0.0, -9.81_f64, 0.0], 1.0 / 60.0_f64, 4);
753        assert_eq!(pipeline.buf.len(), 0);
754        assert!((pipeline.time).abs() < 1e-15);
755    }
756    #[test]
757    fn test_deformable_pipeline_add_particle() {
758        let mut pipeline = DeformableGpuPipeline::new([0.0; 3], 0.016_f64, 2);
759        let idx = pipeline.add_particle([0.0, 5.0_f64, 0.0], 1.0_f64);
760        assert_eq!(idx, 0);
761        assert_eq!(pipeline.buf.len(), 1);
762    }
763    #[test]
764    fn test_deformable_pipeline_step_advances_time() {
765        let mut pipeline = DeformableGpuPipeline::new([0.0, -9.81_f64, 0.0], 0.016_f64, 2);
766        pipeline.add_particle([0.0, 10.0_f64, 0.0], 1.0_f64);
767        pipeline.step();
768        assert!(pipeline.time > 0.0_f64);
769    }
770    #[test]
771    fn test_deformable_pipeline_gravity_falls() {
772        let mut pipeline = DeformableGpuPipeline::new([0.0, -9.81_f64, 0.0], 0.1_f64, 1);
773        pipeline.collision_enabled = false;
774        pipeline.xpbd_enabled = false;
775        pipeline.fem_enabled = false;
776        pipeline.add_particle([0.0, 10.0_f64, 0.0], 1.0_f64);
777        pipeline.step();
778        assert!(
779            pipeline.buf.positions[0][1] < 10.0_f64,
780            "y = {}",
781            pipeline.buf.positions[0][1]
782        );
783    }
784    #[test]
785    fn test_deformable_pipeline_kinetic_energy() {
786        let mut pipeline = DeformableGpuPipeline::new([0.0; 3], 0.016_f64, 1);
787        pipeline.add_particle([0.0; 3], 1.0_f64);
788        assert!((pipeline.kinetic_energy()).abs() < 1e-15);
789    }
790    #[test]
791    fn test_deformable_pipeline_centre_of_mass() {
792        let mut pipeline = DeformableGpuPipeline::new([0.0; 3], 0.016_f64, 1);
793        pipeline.add_particle([0.0; 3], 1.0_f64);
794        pipeline.add_particle([2.0, 0.0, 0.0], 1.0_f64);
795        let com = pipeline.centre_of_mass();
796        assert!((com[0] - 1.0_f64).abs() < 1e-10, "com.x = {}", com[0]);
797    }
798    #[test]
799    fn test_deformable_pipeline_empty_com() {
800        let pipeline = DeformableGpuPipeline::new([0.0; 3], 0.016_f64, 1);
801        let com = pipeline.centre_of_mass();
802        assert_eq!(com, [0.0; 3]);
803    }
804    #[test]
805    fn test_xpbd_solver_with_buffer_distance() {
806        let mut solver = XPBDGpuSolver::new(30);
807        solver.add_distance(0, 1, 1.0_f64, 0.0_f64);
808        let mut buf = DeformableGpuBuffer::new();
809        buf.push([0.0; 3], [0.0; 3], 0.0_f64);
810        buf.push([5.0, 0.0, 0.0], [0.0; 3], 1.0_f64);
811        solver.solve(&mut buf, 0.01_f64);
812        let dist = length3(sub3(buf.positions[1], buf.positions[0]));
813        assert!((dist - 1.0_f64).abs() < 0.01_f64, "dist = {dist}");
814    }
815    #[test]
816    fn test_gpu_buffer_fixed_particle_doesnt_move() {
817        let mut buf = DeformableGpuBuffer::new();
818        buf.push([0.0; 3], [0.0; 3], 0.0_f64);
819        buf.forces[0] = [100.0, 100.0, 100.0];
820        buf.integrate(0.1_f64);
821        assert_eq!(buf.positions[0], [0.0; 3]);
822    }
823}