oxiphysics_softbody/
lib.rs

1// Copyright 2026 COOLJAPAN OU (Team KitaSan)
2// SPDX-License-Identifier: Apache-2.0
3
4//! Soft body simulation for the OxiPhysics engine.
5//!
6//! Provides PBD / XPBD based soft-body simulation including:
7//!
8//! - Distance, bending, volume, and collision constraints
9//! - Cloth simulation with wind forces
10//! - Rope / hair chains
11//! - Corotational FEM soft bodies
12#![warn(missing_docs)]
13
14mod error;
15pub use error::{Error, Result};
16
17pub mod aerodynamics;
18pub mod cloth;
19pub mod constraint;
20pub mod crack_propagation;
21pub mod fem_soft;
22pub mod fracture;
23pub mod fracture_dynamics;
24pub mod inflatable;
25pub mod particle;
26pub mod pbd_system;
27pub mod rope;
28pub mod shape_matching;
29pub mod solver;
30pub mod volume;
31pub mod xpbd;
32
33pub use aerodynamics::AerodynamicsModel;
34pub use cloth::*;
35pub use constraint::{
36    BendingConstraint, CollisionConstraint, DistanceConstraint, SoftConstraint, VolumeConstraint,
37};
38pub use fem_soft::{CorotationalElement, FemSoftBody};
39pub use fracture_dynamics::*;
40pub use inflatable::*;
41pub use particle::{SoftBody, SoftParticle};
42pub use rope::{HairStrand, HairSystem, Rope, RopeSolver, XpbdRope};
43pub use shape_matching::ShapeMatching;
44pub use solver::XpbdSolver;
45
46/// Trait for soft body solvers.
47pub trait SoftBodySolver {
48    /// Initialize this component.
49    fn init(&mut self);
50}
51
52#[cfg(test)]
53mod tests {
54    use super::*;
55    use oxiphysics_core::math::Mat3;
56    use oxiphysics_core::math::{Real, Vec3};
57
58    const EPS: Real = 1e-6;
59
60    // 1. Distance constraint maintains rest length.
61    #[test]
62    fn test_distance_constraint_maintains_rest_length() {
63        let mut particles = vec![
64            SoftParticle::new(Vec3::new(0.0, 0.0, 0.0), 1.0),
65            SoftParticle::new(Vec3::new(2.0, 0.0, 0.0), 1.0),
66        ];
67        let rest = 1.0;
68        let mut dc = DistanceConstraint::new(0, 1, rest, 0.0);
69        // Project many times to converge.
70        for _ in 0..100 {
71            dc.reset_lambda();
72            dc.project(&mut particles, 1.0 / 60.0);
73        }
74        let dist = (particles[0].position - particles[1].position).norm();
75        assert!(
76            (dist - rest).abs() < 0.01,
77            "Distance {dist} should be close to rest length {rest}"
78        );
79    }
80
81    // 2. Two particles connected by distance constraint converge.
82    #[test]
83    fn test_two_particles_converge() {
84        let mut particles = vec![
85            SoftParticle::new(Vec3::new(0.0, 0.0, 0.0), 1.0),
86            SoftParticle::new(Vec3::new(5.0, 0.0, 0.0), 1.0),
87        ];
88        let rest = 1.0;
89        let mut dc = DistanceConstraint::new(0, 1, rest, 0.0);
90        for _ in 0..200 {
91            dc.reset_lambda();
92            dc.project(&mut particles, 1.0 / 60.0);
93        }
94        let dist = (particles[0].position - particles[1].position).norm();
95        assert!(
96            (dist - rest).abs() < 0.05,
97            "Particles should converge to rest length"
98        );
99    }
100
101    // 3. Volume constraint preserves volume.
102    #[test]
103    fn test_volume_constraint_preserves_volume() {
104        let mut particles = vec![
105            SoftParticle::new(Vec3::new(0.0, 0.0, 0.0), 1.0),
106            SoftParticle::new(Vec3::new(1.0, 0.0, 0.0), 1.0),
107            SoftParticle::new(Vec3::new(0.0, 1.0, 0.0), 1.0),
108            SoftParticle::new(Vec3::new(0.0, 0.0, 1.0), 1.0),
109        ];
110        let rest_vol = VolumeConstraint::compute_tet_volume(
111            &particles[0].position,
112            &particles[1].position,
113            &particles[2].position,
114            &particles[3].position,
115        );
116        // Perturb vertex 3.
117        particles[3].position = Vec3::new(0.0, 0.0, 1.5);
118        let mut vc = VolumeConstraint::new([0, 1, 2, 3], rest_vol, 0.0);
119        for _ in 0..100 {
120            vc.reset_lambda();
121            vc.project(&mut particles, 1.0 / 60.0);
122        }
123        let new_vol = VolumeConstraint::compute_tet_volume(
124            &particles[0].position,
125            &particles[1].position,
126            &particles[2].position,
127            &particles[3].position,
128        );
129        assert!(
130            (new_vol - rest_vol).abs() < 0.01,
131            "Volume {new_vol} should be close to rest volume {rest_vol}"
132        );
133    }
134
135    // 4. XPBD solver step doesn't crash.
136    #[test]
137    fn test_xpbd_solver_runs() {
138        let mut body = SoftBody::from_particles(vec![
139            SoftParticle::new(Vec3::new(0.0, 0.0, 0.0), 1.0),
140            SoftParticle::new(Vec3::new(1.0, 0.0, 0.0), 1.0),
141        ]);
142        body.apply_force(&Vec3::new(0.0, -9.81, 0.0));
143        let mut constraints: Vec<Box<dyn SoftConstraint>> = vec![Box::new(
144            DistanceConstraint::from_particles(0, 1, &body.particles, 0.0),
145        )];
146        let mut solver = XpbdSolver::new(5);
147        solver.solve(&mut body, &mut constraints, 1.0 / 60.0);
148        // Just verify no panic and positions have changed.
149        let y0 = body.particles[0].position.y;
150        assert!(y0 < 0.0, "Particle should have fallen under gravity");
151    }
152
153    // 5. Cloth mesh creation has correct topology.
154    #[test]
155    fn test_cloth_mesh_topology() {
156        let cloth = XpbdClothMesh::new(4, 3, 3.0, 2.0, 1.0, 0.001);
157        assert_eq!(cloth.num_particles(), 12); // 4 * 3
158        // Triangles: (nx-1)*(ny-1)*2 = 3*2*2 = 12
159        assert_eq!(cloth.num_triangles(), 12);
160        assert!(!cloth.distance_constraints.is_empty());
161        assert!(!cloth.bending_constraints.is_empty());
162    }
163
164    // 6. Wind force applies in correct direction.
165    #[test]
166    fn test_wind_force_direction() {
167        // A triangle lying in the XZ plane, normal pointing up (+Y).
168        let mut particles = vec![
169            SoftParticle::new(Vec3::new(0.0, 0.0, 0.0), 1.0),
170            SoftParticle::new(Vec3::new(1.0, 0.0, 0.0), 1.0),
171            SoftParticle::new(Vec3::new(0.0, 0.0, 1.0), 1.0),
172        ];
173        let tris = vec![[0, 1, 2]];
174        let wind = Vec3::new(0.0, -1.0, 0.0); // Wind blowing downward.
175        apply_wind(&mut particles, &wind, &tris);
176        // The normal is (0, -1, 0) (depends on winding).  Either way, force
177        // should have a Y component.
178        let total_fy: Real = particles.iter().map(|p| p.external_force.y).sum();
179        assert!(
180            total_fy.abs() > EPS,
181            "Wind should produce a Y-direction force, got {total_fy}"
182        );
183    }
184
185    // 7. Rope creation has correct number of particles/constraints.
186    #[test]
187    fn test_rope_creation() {
188        let rope = XpbdRope::new(
189            Vec3::new(0.0, 10.0, 0.0),
190            Vec3::new(0.0, 0.0, 0.0),
191            10,
192            0.5,
193            0.0,
194        );
195        assert_eq!(rope.num_particles(), 11);
196        assert_eq!(rope.num_segments(), 10);
197        // First particle should be static.
198        assert!(rope.body.particles[0].is_static());
199    }
200
201    // 8. Corotational element identity gives zero force.
202    #[test]
203    fn test_corotational_identity_zero_force() {
204        let particles = vec![
205            SoftParticle::new(Vec3::new(0.0, 0.0, 0.0), 1.0),
206            SoftParticle::new(Vec3::new(1.0, 0.0, 0.0), 1.0),
207            SoftParticle::new(Vec3::new(0.0, 1.0, 0.0), 1.0),
208            SoftParticle::new(Vec3::new(0.0, 0.0, 1.0), 1.0),
209        ];
210        let elem = CorotationalElement::new([0, 1, 2, 3], &particles, 1000.0, 0.3);
211        let forces = elem.compute_forces(&particles);
212        for (k, f) in forces.iter().enumerate() {
213            assert!(
214                f.norm() < 1e-8,
215                "Force on vertex {k} should be zero at rest, got {f:?}"
216            );
217        }
218    }
219
220    // 9. Bending constraint projection.
221    #[test]
222    fn test_bending_constraint_projection() {
223        // Four particles: shared edge (0,1), wing vertices (2,3).
224        let mut particles = vec![
225            SoftParticle::new(Vec3::new(0.0, 0.0, 0.0), 1.0),
226            SoftParticle::new(Vec3::new(1.0, 0.0, 0.0), 1.0),
227            SoftParticle::new(Vec3::new(0.5, 1.0, 0.0), 1.0),
228            SoftParticle::new(Vec3::new(0.5, -1.0, 0.0), 1.0),
229        ];
230        let mut bc = BendingConstraint::from_particles([0, 1, 2, 3], &particles, 0.0);
231        // Perturb wing vertex.
232        particles[2].position = Vec3::new(0.5, 1.0, 0.5);
233        let old_pos = particles[2].position;
234        for _ in 0..50 {
235            bc.reset_lambda();
236            bc.project(&mut particles, 1.0 / 60.0);
237        }
238        // The wing vertex should have moved.
239        let moved = (particles[2].position - old_pos).norm();
240        assert!(
241            moved > 1e-8,
242            "Bending constraint should move perturbed wing vertex"
243        );
244    }
245
246    // 10. Particle free-fall under gravity.
247    #[test]
248    fn test_particle_free_fall() {
249        let mut body =
250            SoftBody::from_particles(vec![SoftParticle::new(Vec3::new(0.0, 10.0, 0.0), 1.0)]);
251        body.apply_force(&Vec3::new(0.0, -9.81, 0.0));
252        let mut constraints: Vec<Box<dyn SoftConstraint>> = Vec::new();
253        let mut solver = XpbdSolver::new(1);
254        let dt = 1.0 / 60.0;
255        for _ in 0..60 {
256            solver.solve(&mut body, &mut constraints, dt);
257        }
258        // After ~1 second of gravity, y should be significantly lower than 10.
259        let y = body.particles[0].position.y;
260        assert!(
261            y < 9.0,
262            "Particle should have fallen under gravity, y = {y}"
263        );
264    }
265
266    // 11. Collision constraint prevents penetration.
267    #[test]
268    fn test_collision_constraint_plane() {
269        let mut particles = vec![SoftParticle::new(Vec3::new(0.0, -1.0, 0.0), 1.0)];
270        let mut cc =
271            CollisionConstraint::new(0, Vec3::new(0.0, 1.0, 0.0), Vec3::new(0.0, 0.0, 0.0));
272        cc.project(&mut particles, 1.0 / 60.0);
273        assert!(
274            particles[0].position.y >= -EPS,
275            "Particle should be pushed above the plane"
276        );
277    }
278
279    // 12. FEM soft body step runs.
280    #[test]
281    fn test_fem_soft_body_step() {
282        let particles = vec![
283            SoftParticle::new_static(Vec3::new(0.0, 0.0, 0.0)),
284            SoftParticle::new(Vec3::new(1.0, 0.0, 0.0), 1.0),
285            SoftParticle::new(Vec3::new(0.0, 1.0, 0.0), 1.0),
286            SoftParticle::new(Vec3::new(0.0, 0.0, 1.0), 1.0),
287        ];
288        let elem = CorotationalElement::new([0, 1, 2, 3], &particles, 1000.0, 0.3);
289        let mut fem = FemSoftBody::new(particles, vec![elem], 0.01);
290        let gravity = Vec3::new(0.0, -9.81, 0.0);
291        for _ in 0..10 {
292            fem.step(1.0 / 60.0, &gravity);
293        }
294        // Dynamic particles should have moved.
295        assert!(
296            fem.particles[1].position != Vec3::new(1.0, 0.0, 0.0),
297            "Dynamic particle should have moved"
298        );
299    }
300
301    // 13. Corotational polar decomposition extracts pure rotation correctly.
302    //
303    // When the deformed tetrahedron is a pure Rz(45°) rotation of the rest
304    // configuration, the deformation gradient F = R (no stretch).  The
305    // corotational model should therefore produce zero elastic forces because
306    // R^T * F - I = R^T * R - I = 0.
307    #[test]
308    fn test_corotational_polar_decomp() {
309        use oxiphysics_core::math::Mat3;
310
311        // Rest tetrahedron.
312        let rest_particles = vec![
313            SoftParticle::new(Vec3::new(0.0, 0.0, 0.0), 1.0),
314            SoftParticle::new(Vec3::new(1.0, 0.0, 0.0), 1.0),
315            SoftParticle::new(Vec3::new(0.0, 1.0, 0.0), 1.0),
316            SoftParticle::new(Vec3::new(0.0, 0.0, 1.0), 1.0),
317        ];
318        let elem = CorotationalElement::new([0, 1, 2, 3], &rest_particles, 1000.0, 0.3);
319
320        // Apply a pure Rz(45°) rotation to all vertices.
321        let angle: Real = std::f64::consts::FRAC_PI_4; // 45 degrees
322        let cos_a = angle.cos();
323        let sin_a = angle.sin();
324        let rz = Mat3::new(cos_a, -sin_a, 0.0, sin_a, cos_a, 0.0, 0.0, 0.0, 1.0);
325        let rotated_particles: Vec<SoftParticle> = rest_particles
326            .iter()
327            .map(|p| SoftParticle::new(rz * p.position, 1.0))
328            .collect();
329
330        // Forces on a purely-rotated (unstretched) element must be zero.
331        let forces = elem.compute_forces(&rotated_particles);
332        for (k, f) in forces.iter().enumerate() {
333            assert!(
334                f.norm() < 1e-8,
335                "Force on vertex {k} should be zero for pure rotation, got {f:?}"
336            );
337        }
338    }
339
340    // 14. Cloth sags under gravity (center vertex moves downward).
341    //
342    // A 3x3 grid cloth is pinned at its four corners.  After 100 XPBD steps
343    // with gravity applied, the center particle (index 4) must have a negative
344    // Y position (sag below the initial flat XZ plane).
345    #[test]
346    fn test_cloth_gravity_sag() {
347        let mut cloth = XpbdClothMesh::new(3, 3, 2.0, 2.0, 1.0, 1e-4);
348        // Pin the four corners.
349        cloth.pin(0, 0);
350        cloth.pin(2, 0);
351        cloth.pin(0, 2);
352        cloth.pin(2, 2);
353
354        let gravity = Vec3::new(0.0, -9.81, 0.0);
355        let dt = 1.0 / 60.0;
356        let mut solver = XpbdSolver::new(10);
357
358        // Collect all constraints into a single boxed list.
359        let mut constraints: Vec<Box<dyn SoftConstraint>> = cloth
360            .distance_constraints
361            .iter()
362            .cloned()
363            .map(|c| Box::new(c) as Box<dyn SoftConstraint>)
364            .chain(
365                cloth
366                    .bending_constraints
367                    .iter()
368                    .cloned()
369                    .map(|c| Box::new(c) as Box<dyn SoftConstraint>),
370            )
371            .collect();
372
373        for _ in 0..100 {
374            cloth.body.apply_force(&gravity);
375            solver.solve(&mut cloth.body, &mut constraints, dt);
376            cloth.body.clear_forces();
377        }
378
379        // Center particle is index nx/2 * nx + nx/2 = 1 * 3 + 1 = 4 for a 3x3 grid.
380        let center_y = cloth.body.particles[4].position.y;
381        assert!(
382            center_y < -0.01,
383            "Center particle should sag below y=0 under gravity, got y={center_y}"
384        );
385    }
386
387    // -----------------------------------------------------------------------
388    // Aerodynamics tests
389    // -----------------------------------------------------------------------
390
391    // A1. Zero wind velocity → zero forces on every vertex.
392    #[test]
393    fn test_aero_zero_wind() {
394        let aero = AerodynamicsModel::new(1.225, 1.0, 0.5);
395        let zero = Vec3::zeros();
396        let (f0, f1, f2) = aero.triangle_force(
397            Vec3::new(0.0, 0.0, 0.0),
398            Vec3::new(1.0, 0.0, 0.0),
399            Vec3::new(0.0, 0.0, 1.0),
400            zero,
401            zero,
402            zero,
403            zero, // zero wind
404        );
405        assert!(
406            f0.norm() < 1e-10,
407            "f0 should be zero with zero wind, got {:?}",
408            f0
409        );
410        assert!(
411            f1.norm() < 1e-10,
412            "f1 should be zero with zero wind, got {:?}",
413            f1
414        );
415        assert!(
416            f2.norm() < 1e-10,
417            "f2 should be zero with zero wind, got {:?}",
418            f2
419        );
420    }
421
422    // A2. Wind perpendicular to triangle (head-on) → purely drag force.
423    //
424    // Triangle in XZ plane (normal along +Y).  Wind blows in +Y direction.
425    // The lift term vanishes because v_rel × n̂ = 0 when v_rel ∥ n̂.
426    #[test]
427    fn test_aero_head_on_wind() {
428        let cd = 1.5;
429        let cl = 0.8;
430        let aero = AerodynamicsModel::new(1.225, cd, cl);
431        let zero = Vec3::zeros();
432        let wind = Vec3::new(0.0, 10.0, 0.0); // perpendicular to XZ plane
433        let (f0, f1, f2) = aero.triangle_force(
434            Vec3::new(0.0, 0.0, 0.0),
435            Vec3::new(1.0, 0.0, 0.0),
436            Vec3::new(0.0, 0.0, 1.0),
437            zero,
438            zero,
439            zero,
440            wind,
441        );
442        let f_total = f0 + f1 + f2;
443        // Force must be entirely in the Y direction (within tolerance).
444        assert!(
445            f_total.x.abs() < 1e-10 && f_total.z.abs() < 1e-10,
446            "Head-on wind should produce only Y force, got {:?}",
447            f_total
448        );
449        assert!(
450            f_total.y > 1e-10,
451            "Head-on wind should produce positive Y drag, got {:?}",
452            f_total
453        );
454    }
455
456    // A3. Doubling rho_air doubles the total aerodynamic force.
457    #[test]
458    fn test_aero_force_scales_with_rho() {
459        let aero1 = AerodynamicsModel::new(1.0, 1.0, 0.5);
460        let aero2 = AerodynamicsModel::new(2.0, 1.0, 0.5);
461        let zero = Vec3::zeros();
462        let wind = Vec3::new(5.0, 3.0, 1.0);
463        let v0 = Vec3::new(0.0, 0.0, 0.0);
464        let v1 = Vec3::new(1.0, 0.0, 0.0);
465        let v2 = Vec3::new(0.0, 0.0, 1.0);
466
467        let (a0, a1, a2) = aero1.triangle_force(v0, v1, v2, zero, zero, zero, wind);
468        let (b0, b1, b2) = aero2.triangle_force(v0, v1, v2, zero, zero, zero, wind);
469
470        let fa = (a0 + a1 + a2).norm();
471        let fb = (b0 + b1 + b2).norm();
472
473        assert!(
474            (fb - 2.0 * fa).abs() < 1e-10,
475            "Doubling rho should double force: fa={fa}, fb={fb}"
476        );
477    }
478
479    // -----------------------------------------------------------------------
480    // Shape matching tests
481    // -----------------------------------------------------------------------
482
483    // S1. Pure rotation: shape-matching goal positions = rotated rest positions.
484    //
485    // If the current positions are a pure rotation R of the rest positions
486    // (with stiffness=1), the goals must equal the current positions.
487    #[test]
488    fn test_shape_matching_rigid_body() {
489        let rest = vec![
490            Vec3::new(1.0, 0.0, 0.0),
491            Vec3::new(-1.0, 0.0, 0.0),
492            Vec3::new(0.0, 1.0, 0.0),
493            Vec3::new(0.0, -1.0, 0.0),
494        ];
495        let masses = vec![1.0f64; 4];
496        let sm = ShapeMatching::new(&rest, &masses, 1.0);
497
498        // Apply a 90° rotation about Z.
499        let rz90 = Mat3::new(0.0, -1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0);
500        let current: Vec<Vec3> = rest.iter().map(|p| rz90 * p).collect();
501        let goals = sm.goal_positions(&current);
502
503        for (i, (goal, cur)) in goals.iter().zip(current.iter()).enumerate() {
504            assert!(
505                (goal - cur).norm() < 1e-6,
506                "Particle {i}: goal {:?} should equal current {:?}",
507                goal,
508                cur
509            );
510        }
511    }
512
513    // S2. stiffness=0 → goals equal current positions (no constraint).
514    #[test]
515    fn test_shape_matching_stiffness_zero() {
516        let rest = vec![
517            Vec3::new(1.0, 0.0, 0.0),
518            Vec3::new(-1.0, 0.0, 0.0),
519            Vec3::new(0.0, 2.0, 0.0),
520        ];
521        let masses = vec![1.0f64; 3];
522        let sm = ShapeMatching::new(&rest, &masses, 0.0); // stiffness = 0
523
524        // Arbitrary deformed positions.
525        let current = vec![
526            Vec3::new(3.0, 1.0, -2.0),
527            Vec3::new(-2.0, 0.5, 1.0),
528            Vec3::new(0.5, 4.0, 0.3),
529        ];
530        let goals = sm.goal_positions(&current);
531
532        for (i, (goal, cur)) in goals.iter().zip(current.iter()).enumerate() {
533            assert!(
534                (goal - cur).norm() < 1e-10,
535                "With stiffness=0 goal {i} {:?} should equal current {:?}",
536                goal,
537                cur
538            );
539        }
540    }
541
542    // S3. Identity: undeformed positions → goals equal current positions.
543    #[test]
544    fn test_shape_matching_identity() {
545        let rest = vec![
546            Vec3::new(0.0, 0.0, 0.0),
547            Vec3::new(1.0, 0.0, 0.0),
548            Vec3::new(0.0, 1.0, 0.0),
549            Vec3::new(0.0, 0.0, 1.0),
550        ];
551        let masses = vec![1.0f64; 4];
552        let sm = ShapeMatching::new(&rest, &masses, 1.0);
553
554        // Current positions = rest positions (no deformation).
555        let goals = sm.goal_positions(&rest);
556
557        for (i, (goal, cur)) in goals.iter().zip(rest.iter()).enumerate() {
558            assert!(
559                (goal - cur).norm() < 1e-6,
560                "Particle {i}: goal {:?} should equal rest {:?}",
561                goal,
562                cur
563            );
564        }
565    }
566
567    // 15. Distance constraint pushes two particles that are too close to rest length.
568    //
569    // Place two particles at the same position (distance = 0).  After many
570    // XPBD iterations the distance must converge to the rest length.
571    #[test]
572    fn test_distance_constraint_correction() {
573        let mut particles = vec![
574            SoftParticle::new(Vec3::new(0.0, 0.0, 0.0), 1.0),
575            SoftParticle::new(Vec3::new(0.0, 0.0, 0.0), 1.0), // same position
576        ];
577        let rest = 1.5;
578        // To avoid the degenerate (dist < 1e-12) guard, give them a tiny offset.
579        particles[1].position = Vec3::new(1e-6, 0.0, 0.0);
580
581        let mut dc = DistanceConstraint::new(0, 1, rest, 0.0);
582        for _ in 0..200 {
583            dc.reset_lambda();
584            dc.project(&mut particles, 1.0 / 60.0);
585        }
586        let dist = (particles[0].position - particles[1].position).norm();
587        assert!(
588            (dist - rest).abs() < 0.1,
589            "Particles too close should be pushed to rest length {rest}, got dist={dist}"
590        );
591    }
592
593    // 16. Volume constraint keeps tetrahedron volume close to rest after perturbation.
594    //
595    // Form a unit tetrahedron, record its volume, then compress one vertex
596    // inward.  After running the volume constraint the volume should return
597    // to within 5 % of the original.
598    #[test]
599    fn test_volume_constraint_preserve() {
600        let mut particles = vec![
601            SoftParticle::new(Vec3::new(0.0, 0.0, 0.0), 1.0),
602            SoftParticle::new(Vec3::new(2.0, 0.0, 0.0), 1.0),
603            SoftParticle::new(Vec3::new(0.0, 2.0, 0.0), 1.0),
604            SoftParticle::new(Vec3::new(0.0, 0.0, 2.0), 1.0),
605        ];
606        let rest_vol = VolumeConstraint::compute_tet_volume(
607            &particles[0].position,
608            &particles[1].position,
609            &particles[2].position,
610            &particles[3].position,
611        );
612
613        // Compress: move vertex 3 so volume becomes roughly half.
614        particles[3].position = Vec3::new(0.0, 0.0, 0.5);
615
616        let mut vc = VolumeConstraint::new([0, 1, 2, 3], rest_vol, 0.0);
617        for _ in 0..200 {
618            vc.reset_lambda();
619            vc.project(&mut particles, 1.0 / 60.0);
620        }
621
622        let final_vol = VolumeConstraint::compute_tet_volume(
623            &particles[0].position,
624            &particles[1].position,
625            &particles[2].position,
626            &particles[3].position,
627        );
628        let relative_error = ((final_vol - rest_vol) / rest_vol).abs();
629        assert!(
630            relative_error < 0.05,
631            "Volume after constraint should be close to rest ({rest_vol}), got {final_vol} (error={relative_error:.3})"
632        );
633    }
634}
635
636pub mod active_matter;
637pub mod active_origami;
638pub mod bio_softbody;
639pub mod bioinspired;
640pub mod biomech_simulation;
641pub mod cable_nets;
642pub mod cloth_advanced;
643pub mod cloth_simulation;
644pub mod collision_response;
645pub mod constraint_visualization;
646pub mod cosserat_rods;
647pub mod cosserat_softbody;
648pub mod elastic_wave;
649pub mod electroactive_softbody;
650pub mod fluid_film_softbody;
651pub mod food_physics;
652pub mod granular_softbody;
653pub mod growth_mechanics;
654pub mod hair;
655pub mod hair_fur;
656pub mod hair_sim;
657pub mod hair_softbody;
658pub mod haptic_softbody;
659pub mod liquid_sim;
660pub mod material_point;
661pub mod membrane_biophysics;
662pub mod membrane_softbody;
663pub mod metamaterial_softbody;
664pub mod metamorphic_mesh;
665pub mod morphogenesis_softbody;
666pub mod mud_snow;
667pub mod muscle_sim;
668pub mod muscle_simulation;
669pub mod muscle_softbody;
670pub mod muscle_tendon;
671pub mod neural_deform;
672pub mod neural_soft;
673pub mod neural_softbody;
674pub mod numerical_softbody;
675pub mod origami_folding;
676pub mod origami_mech;
677pub mod origami_mechanics;
678pub mod origami_softbody;
679pub mod pbd_constraints;
680pub mod peridynamics;
681pub mod position_based_fluids;
682pub mod pressure_volume;
683pub mod rigid_soft_coupling;
684pub mod rod_mechanics;
685pub mod rope_simulation;
686pub mod sand_sim;
687pub mod smart_materials;
688pub mod soft_gripper;
689pub mod surgery_simulation;
690pub mod surgical_simulation;
691pub mod tearing;
692pub mod tendon_softbody;
693pub mod textile_simulation;
694pub mod tissue_sim;
695pub mod topology_opt;
696pub mod topology_optimization;
697pub mod underwater_softbody;
698pub mod wrinkling;
oxiphysics_softbody/lib.rs

oxiphysics_softbody/
lib.rs
