From 70afc5a7d871919776a64782e8b93404e6b0defd Mon Sep 17 00:00:00 2001 From: eug-vs Date: Sun, 15 Dec 2024 13:17:43 +0100 Subject: feat!: add raylib rendering --- .gitignore | 2 - Cargo.lock | 201 ---------- Cargo.toml | 7 - README.md | 4 - physics/.gitignore | 2 + physics/Cargo.lock | 201 ++++++++++ physics/Cargo.toml | 7 + physics/README.md | 4 + physics/src/algebra/mod.rs | 9 + physics/src/algebra/subspace.rs | 70 ++++ physics/src/constraint/anchor.rs | 43 +++ physics/src/constraint/beam.rs | 47 +++ physics/src/constraint/mod.rs | 140 +++++++ physics/src/constraint/slider.rs | 47 +++ physics/src/force/drag.rs | 15 + physics/src/force/gravity.rs | 15 + physics/src/force/mod.rs | 24 ++ physics/src/force/spring.rs | 28 ++ physics/src/lib.rs | 9 + physics/src/particle_system.rs | 53 +++ physics/src/ppm.rs | 67 ++++ physics/src/renderer/mod.rs | 68 ++++ physics/src/solver/midpoint.rs | 18 + physics/src/solver/mod.rs | 188 ++++++++++ playground/.gitignore | 2 + playground/Cargo.lock | 783 +++++++++++++++++++++++++++++++++++++++ playground/Cargo.toml | 8 + playground/src/main.rs | 185 +++++++++ src/algebra/mod.rs | 1 - src/algebra/subspace.rs | 69 ---- src/constraint/anchor.rs | 42 --- src/constraint/beam.rs | 46 --- src/constraint/mod.rs | 139 ------- src/constraint/slider.rs | 50 --- src/force/drag.rs | 15 - src/force/gravity.rs | 15 - src/force/mod.rs | 24 -- src/main.rs | 59 --- src/particle_system.rs | 47 --- src/ppm.rs | 41 -- src/solver/midpoint.rs | 18 - src/solver/mod.rs | 187 ---------- 42 files changed, 2033 insertions(+), 967 deletions(-) delete mode 100644 .gitignore delete mode 100644 Cargo.lock delete mode 100644 Cargo.toml delete mode 100644 README.md create mode 100644 physics/.gitignore create mode 100644 physics/Cargo.lock create mode 100644 physics/Cargo.toml create mode 100644 physics/README.md create mode 100644 physics/src/algebra/mod.rs create mode 100644 physics/src/algebra/subspace.rs create mode 100644 physics/src/constraint/anchor.rs create mode 100644 physics/src/constraint/beam.rs create mode 100644 physics/src/constraint/mod.rs create mode 100644 physics/src/constraint/slider.rs create mode 100644 physics/src/force/drag.rs create mode 100644 physics/src/force/gravity.rs create mode 100644 physics/src/force/mod.rs create mode 100644 physics/src/force/spring.rs create mode 100644 physics/src/lib.rs create mode 100644 physics/src/particle_system.rs create mode 100644 physics/src/ppm.rs create mode 100644 physics/src/renderer/mod.rs create mode 100644 physics/src/solver/midpoint.rs create mode 100644 physics/src/solver/mod.rs create mode 100644 playground/.gitignore create mode 100644 playground/Cargo.lock create mode 100644 playground/Cargo.toml create mode 100644 playground/src/main.rs delete mode 100644 src/algebra/mod.rs delete mode 100644 src/algebra/subspace.rs delete mode 100644 src/constraint/anchor.rs delete mode 100644 src/constraint/beam.rs delete mode 100644 src/constraint/mod.rs delete mode 100644 src/constraint/slider.rs delete mode 100644 src/force/drag.rs delete mode 100644 src/force/gravity.rs delete mode 100644 src/force/mod.rs delete mode 100644 src/main.rs delete mode 100644 src/particle_system.rs delete mode 100644 src/ppm.rs delete mode 100644 src/solver/midpoint.rs delete mode 100644 src/solver/mod.rs diff --git a/.gitignore b/.gitignore deleted file mode 100644 index c868ffa..0000000 --- a/.gitignore +++ /dev/null @@ -1,2 +0,0 @@ -/target -out diff --git a/Cargo.lock b/Cargo.lock deleted file mode 100644 index a26bd90..0000000 --- a/Cargo.lock +++ /dev/null @@ -1,201 +0,0 @@ -# This file is automatically @generated by Cargo. -# It is not intended for manual editing. -version = 3 - -[[package]] -name = "approx" -version = "0.5.1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "cab112f0a86d568ea0e627cc1d6be74a1e9cd55214684db5561995f6dad897c6" -dependencies = [ - "num-traits", -] - -[[package]] -name = "autocfg" -version = "1.4.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "ace50bade8e6234aa140d9a2f552bbee1db4d353f69b8217bc503490fc1a9f26" - -[[package]] -name = "bytemuck" -version = "1.20.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "8b37c88a63ffd85d15b406896cc343916d7cf57838a847b3a6f2ca5d39a5695a" - -[[package]] -name = "matrixmultiply" -version = "0.3.9" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "9380b911e3e96d10c1f415da0876389aaf1b56759054eeb0de7df940c456ba1a" -dependencies = [ - "autocfg", - "rawpointer", -] - -[[package]] -name = "nalgebra" -version = "0.33.2" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "26aecdf64b707efd1310e3544d709c5c0ac61c13756046aaaba41be5c4f66a3b" -dependencies = [ - "approx", - "matrixmultiply", - "nalgebra-macros", - "num-complex", - "num-rational", - "num-traits", - "simba", - "typenum", -] - -[[package]] -name = "nalgebra-macros" -version = "0.2.2" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "254a5372af8fc138e36684761d3c0cdb758a4410e938babcff1c860ce14ddbfc" -dependencies = [ - "proc-macro2", - "quote", - "syn", -] - -[[package]] -name = "num-bigint" -version = "0.4.6" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "a5e44f723f1133c9deac646763579fdb3ac745e418f2a7af9cd0c431da1f20b9" -dependencies = [ - "num-integer", - "num-traits", -] - -[[package]] -name = "num-complex" -version = "0.4.6" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "73f88a1307638156682bada9d7604135552957b7818057dcef22705b4d509495" -dependencies = [ - "num-traits", -] - -[[package]] -name = "num-integer" -version = "0.1.46" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "7969661fd2958a5cb096e56c8e1ad0444ac2bbcd0061bd28660485a44879858f" -dependencies = [ - "num-traits", -] - -[[package]] -name = "num-rational" -version = "0.4.2" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "f83d14da390562dca69fc84082e73e548e1ad308d24accdedd2720017cb37824" -dependencies = [ - "num-bigint", - "num-integer", - "num-traits", -] - -[[package]] -name = "num-traits" -version = "0.2.19" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "071dfc062690e90b734c0b2273ce72ad0ffa95f0c74596bc250dcfd960262841" -dependencies = [ - "autocfg", -] - -[[package]] -name = "paste" -version = "1.0.15" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "57c0d7b74b563b49d38dae00a0c37d4d6de9b432382b2892f0574ddcae73fd0a" - -[[package]] -name = "physics-rust" -version = "0.1.0" -dependencies = [ - "nalgebra", -] - -[[package]] -name = "proc-macro2" -version = "1.0.92" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "37d3544b3f2748c54e147655edb5025752e2303145b5aefb3c3ea2c78b973bb0" -dependencies = [ - "unicode-ident", -] - -[[package]] -name = "quote" -version = "1.0.37" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "b5b9d34b8991d19d98081b46eacdd8eb58c6f2b201139f7c5f643cc155a633af" -dependencies = [ - "proc-macro2", -] - -[[package]] -name = "rawpointer" -version = "0.2.1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "60a357793950651c4ed0f3f52338f53b2f809f32d83a07f72909fa13e4c6c1e3" - -[[package]] -name = "safe_arch" -version = "0.7.2" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "c3460605018fdc9612bce72735cba0d27efbcd9904780d44c7e3a9948f96148a" -dependencies = [ - "bytemuck", -] - -[[package]] -name = "simba" -version = "0.9.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "b3a386a501cd104797982c15ae17aafe8b9261315b5d07e3ec803f2ea26be0fa" -dependencies = [ - "approx", - "num-complex", - "num-traits", - "paste", - "wide", -] - -[[package]] -name = "syn" -version = "2.0.90" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "919d3b74a5dd0ccd15aeb8f93e7006bd9e14c295087c9896a110f490752bcf31" -dependencies = [ - "proc-macro2", - "quote", - "unicode-ident", -] - -[[package]] -name = "typenum" -version = "1.17.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "42ff0bf0c66b8238c6f3b578df37d0b7848e55df8577b3f74f92a69acceeb825" - -[[package]] -name = "unicode-ident" -version = "1.0.14" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "adb9e6ca4f869e1180728b7950e35922a7fc6397f7b641499e8f3ef06e50dc83" - -[[package]] -name = "wide" -version = "0.7.30" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "58e6db2670d2be78525979e9a5f9c69d296fd7d670549fe9ebf70f8708cb5019" -dependencies = [ - "bytemuck", - "safe_arch", -] diff --git a/Cargo.toml b/Cargo.toml deleted file mode 100644 index 11c9611..0000000 --- a/Cargo.toml +++ /dev/null @@ -1,7 +0,0 @@ -[package] -name = "physics-rust" -version = "0.1.0" -edition = "2021" - -[dependencies] -nalgebra = "0.33.2" diff --git a/README.md b/README.md deleted file mode 100644 index 59d5c16..0000000 --- a/README.md +++ /dev/null @@ -1,4 +0,0 @@ -# Useful links -## Constrained dynamics - - [A lot of pictures to get general idea and familiarize, but not very understandable](https://danielchappuis.ch/download/ConstraintsDerivationRigidBody3D.pdf) - - [More text and math but very well written and understandable](https://graphics.pixar.com/pbm2001/pdf/notesf.pdf) diff --git a/physics/.gitignore b/physics/.gitignore new file mode 100644 index 0000000..c868ffa --- /dev/null +++ b/physics/.gitignore @@ -0,0 +1,2 @@ +/target +out diff --git a/physics/Cargo.lock b/physics/Cargo.lock new file mode 100644 index 0000000..98fe5af --- /dev/null +++ b/physics/Cargo.lock @@ -0,0 +1,201 @@ +# This file is automatically @generated by Cargo. +# It is not intended for manual editing. +version = 3 + +[[package]] +name = "approx" +version = "0.5.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "cab112f0a86d568ea0e627cc1d6be74a1e9cd55214684db5561995f6dad897c6" +dependencies = [ + "num-traits", +] + +[[package]] +name = "autocfg" +version = "1.4.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ace50bade8e6234aa140d9a2f552bbee1db4d353f69b8217bc503490fc1a9f26" + +[[package]] +name = "bytemuck" +version = "1.20.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8b37c88a63ffd85d15b406896cc343916d7cf57838a847b3a6f2ca5d39a5695a" + +[[package]] +name = "matrixmultiply" +version = "0.3.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9380b911e3e96d10c1f415da0876389aaf1b56759054eeb0de7df940c456ba1a" +dependencies = [ + "autocfg", + "rawpointer", +] + +[[package]] +name = "nalgebra" +version = "0.33.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "26aecdf64b707efd1310e3544d709c5c0ac61c13756046aaaba41be5c4f66a3b" +dependencies = [ + "approx", + "matrixmultiply", + "nalgebra-macros", + "num-complex", + "num-rational", + "num-traits", + "simba", + "typenum", +] + +[[package]] +name = "nalgebra-macros" +version = "0.2.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "254a5372af8fc138e36684761d3c0cdb758a4410e938babcff1c860ce14ddbfc" +dependencies = [ + "proc-macro2", + "quote", + "syn", +] + +[[package]] +name = "num-bigint" +version = "0.4.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a5e44f723f1133c9deac646763579fdb3ac745e418f2a7af9cd0c431da1f20b9" +dependencies = [ + "num-integer", + "num-traits", +] + +[[package]] +name = "num-complex" +version = "0.4.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "73f88a1307638156682bada9d7604135552957b7818057dcef22705b4d509495" +dependencies = [ + "num-traits", +] + +[[package]] +name = "num-integer" +version = "0.1.46" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7969661fd2958a5cb096e56c8e1ad0444ac2bbcd0061bd28660485a44879858f" +dependencies = [ + "num-traits", +] + +[[package]] +name = "num-rational" +version = "0.4.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f83d14da390562dca69fc84082e73e548e1ad308d24accdedd2720017cb37824" +dependencies = [ + "num-bigint", + "num-integer", + "num-traits", +] + +[[package]] +name = "num-traits" +version = "0.2.19" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "071dfc062690e90b734c0b2273ce72ad0ffa95f0c74596bc250dcfd960262841" +dependencies = [ + "autocfg", +] + +[[package]] +name = "paste" +version = "1.0.15" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "57c0d7b74b563b49d38dae00a0c37d4d6de9b432382b2892f0574ddcae73fd0a" + +[[package]] +name = "physics" +version = "0.1.0" +dependencies = [ + "nalgebra", +] + +[[package]] +name = "proc-macro2" +version = "1.0.92" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "37d3544b3f2748c54e147655edb5025752e2303145b5aefb3c3ea2c78b973bb0" +dependencies = [ + "unicode-ident", +] + +[[package]] +name = "quote" +version = "1.0.37" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b5b9d34b8991d19d98081b46eacdd8eb58c6f2b201139f7c5f643cc155a633af" +dependencies = [ + "proc-macro2", +] + +[[package]] +name = "rawpointer" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "60a357793950651c4ed0f3f52338f53b2f809f32d83a07f72909fa13e4c6c1e3" + +[[package]] +name = "safe_arch" +version = "0.7.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c3460605018fdc9612bce72735cba0d27efbcd9904780d44c7e3a9948f96148a" +dependencies = [ + "bytemuck", +] + +[[package]] +name = "simba" +version = "0.9.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b3a386a501cd104797982c15ae17aafe8b9261315b5d07e3ec803f2ea26be0fa" +dependencies = [ + "approx", + "num-complex", + "num-traits", + "paste", + "wide", +] + +[[package]] +name = "syn" +version = "2.0.90" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "919d3b74a5dd0ccd15aeb8f93e7006bd9e14c295087c9896a110f490752bcf31" +dependencies = [ + "proc-macro2", + "quote", + "unicode-ident", +] + +[[package]] +name = "typenum" +version = "1.17.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "42ff0bf0c66b8238c6f3b578df37d0b7848e55df8577b3f74f92a69acceeb825" + +[[package]] +name = "unicode-ident" +version = "1.0.14" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "adb9e6ca4f869e1180728b7950e35922a7fc6397f7b641499e8f3ef06e50dc83" + +[[package]] +name = "wide" +version = "0.7.30" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "58e6db2670d2be78525979e9a5f9c69d296fd7d670549fe9ebf70f8708cb5019" +dependencies = [ + "bytemuck", + "safe_arch", +] diff --git a/physics/Cargo.toml b/physics/Cargo.toml new file mode 100644 index 0000000..e5db63e --- /dev/null +++ b/physics/Cargo.toml @@ -0,0 +1,7 @@ +[package] +name = "physics" +version = "0.1.0" +edition = "2021" + +[dependencies] +nalgebra = "0.33.2" diff --git a/physics/README.md b/physics/README.md new file mode 100644 index 0000000..59d5c16 --- /dev/null +++ b/physics/README.md @@ -0,0 +1,4 @@ +# Useful links +## Constrained dynamics + - [A lot of pictures to get general idea and familiarize, but not very understandable](https://danielchappuis.ch/download/ConstraintsDerivationRigidBody3D.pdf) + - [More text and math but very well written and understandable](https://graphics.pixar.com/pbm2001/pdf/notesf.pdf) diff --git a/physics/src/algebra/mod.rs b/physics/src/algebra/mod.rs new file mode 100644 index 0000000..47efba8 --- /dev/null +++ b/physics/src/algebra/mod.rs @@ -0,0 +1,9 @@ +use nalgebra::{Point as PointBase, SVector}; + +pub const N: usize = 2; +pub type Scalar = f64; + +pub type Vector = SVector; +pub type Point = PointBase; + +pub mod subspace; diff --git a/physics/src/algebra/subspace.rs b/physics/src/algebra/subspace.rs new file mode 100644 index 0000000..9fb944e --- /dev/null +++ b/physics/src/algebra/subspace.rs @@ -0,0 +1,70 @@ +use nalgebra::SMatrix; + +use super::{Scalar, N, Point, Vector}; + + +type ProjectionMatrix = SMatrix; + +pub struct Subspace { + point: Point, + vectors: [Vector; DIM], + projection_matrix: ProjectionMatrix, +} + +pub type Line = Subspace<1>; +pub type Plane = Subspace<2>; + +impl Subspace { + pub fn new(point: Point, mut vectors: [Vector; DIM]) -> Self { + for vector in &mut vectors { + vector.normalize_mut(); + } + + Self { + point, + vectors, + projection_matrix: { + let v = SMatrix::::from_columns(&vectors); + let transpose = v.transpose(); + let inverse = (transpose * v).try_inverse().unwrap(); + v * inverse * transpose + } + } + } + + pub fn project_point(&self, point: Point) -> Point { + self.projection_matrix * (point - self.point.coords) + self.point.coords + } + + pub fn distance(&self, point: Point) -> Scalar { + let projected = self.project_point(point); + (projected - point).norm() + } +} + +#[cfg(test)] +mod tests { + use crate::algebra::{subspace::Line, Point, Vector}; + + #[test] + fn test_projection_onto_line() { + let line = Line::new(Point::new(1.0, 0.0), [Vector::new(3.0, 1.0)]); + let point = Point::new(3.0, 4.0); + + let matrix = line.projection_matrix; + { + let diff = matrix * matrix - matrix; + assert!(diff.norm() < 0.001, "Projection matrix squared should be equal to itself: {}, {}", matrix, matrix* matrix); + } + { + let diff = matrix - matrix.transpose(); + assert!(diff.norm() < 0.001, "Projection matrix transpose should be equal to itself: {}, {}", matrix, matrix.transpose()); + } + + { + let projected = line.project_point(point); + let diff = projected - Point::new(4.0, 1.0); + assert!(diff.norm() < 0.001, "Point projected incorrectly"); + } + } +} diff --git a/physics/src/constraint/anchor.rs b/physics/src/constraint/anchor.rs new file mode 100644 index 0000000..c90cfb2 --- /dev/null +++ b/physics/src/constraint/anchor.rs @@ -0,0 +1,43 @@ +use nalgebra::{DVector, RowDVector}; + +use crate::particle_system::ParticleSystem; +use crate::algebra::{Point, Scalar, N}; + +use super::Constraint; + +pub struct AnchorConstraint { + pub particle_id: usize, + pub anchor: Point, + + jacobian: RowDVector, +} + +impl ParticleSystem { + pub fn add_anchor_constraint(&mut self, particle_id: usize) { + let anchor = self.particles[particle_id].position; + self.constraints.push(Box::new(AnchorConstraint { + particle_id, + anchor, + jacobian: RowDVector::zeros(self.particles.len() * N), + })); + } +} + +impl Constraint for AnchorConstraint { + fn get_particles(&self) -> Vec { + vec![self.particle_id] + } + + fn c(&self, q: &DVector) -> Scalar { + let position = q.fixed_rows(self.particle_id * N); + (position - self.anchor.coords).norm() + } + + fn set_jacobian(&mut self, jacobian: RowDVector) { + self.jacobian = jacobian + } + + fn jacobian_prev(&self) -> RowDVector { + self.jacobian.clone() + } +} diff --git a/physics/src/constraint/beam.rs b/physics/src/constraint/beam.rs new file mode 100644 index 0000000..14b1c1f --- /dev/null +++ b/physics/src/constraint/beam.rs @@ -0,0 +1,47 @@ +use nalgebra::{DVector, RowDVector}; + +use crate::particle_system::ParticleSystem; +use crate::algebra::{Scalar, N}; + +use super::Constraint; + +pub struct BeamConstraint { + pub particle_ids: [usize; 2], + pub length: Scalar, + + jacobian: RowDVector, +} + +impl ParticleSystem { + pub fn add_beam_constraint(&mut self, particle_ids: [usize; 2]) { + let a = &self.particles[particle_ids[0]]; + let b = &self.particles[particle_ids[1]]; + + self.constraints.push(Box::new(BeamConstraint { + particle_ids, + length: (a.position - b.position).norm(), + jacobian: RowDVector::zeros(self.particles.len() * N), + })); + } +} + +impl Constraint for BeamConstraint { + fn get_particles(&self) -> Vec { + Vec::from(self.particle_ids) + } + + fn c(&self, q: &DVector) -> Scalar { + let a = q.fixed_rows::(self.particle_ids[0] * N); + let b = q.fixed_rows::(self.particle_ids[1] * N); + + (a - b).norm() - self.length + } + + fn set_jacobian(&mut self, jacobian: RowDVector) { + self.jacobian = jacobian + } + + fn jacobian_prev(&self) -> RowDVector { + self.jacobian.clone() + } +} diff --git a/physics/src/constraint/mod.rs b/physics/src/constraint/mod.rs new file mode 100644 index 0000000..b2958ec --- /dev/null +++ b/physics/src/constraint/mod.rs @@ -0,0 +1,140 @@ +use nalgebra::{DMatrix, DVector, RowDVector}; + +use crate::particle_system::ParticleSystem; +use crate::algebra::{Scalar, Vector, N}; +pub mod anchor; +pub mod slider; +pub mod beam; + +const SPRING_CONSTANT: Scalar = 0.65; +const DAMPING_CONSTANT: Scalar = 0.85; + +/// SIZE is always 3xN, but this operation can't be done at compile time yet: +/// "generic parameters may not be used in const operations" +pub trait Constraint { + fn get_particles(&self) -> Vec; + + fn c(&self, q: &DVector) -> Scalar; + + fn jacobian_prev(&self) -> RowDVector; + fn set_jacobian(&mut self, jacobian: RowDVector); + + fn partial_derivative(&self, q: &DVector) -> RowDVector { + let dq = 0.001; + let c_original = self.c(&q); + let mut result = RowDVector::zeros(q.len()); + // The only non-zero components of derivative vector are + // the particles that this constraint affects + for particle_id in self.get_particles() { + for i in 0..N { + let index = particle_id * N + i; + let mut q_partial = q.clone(); + q_partial[index] += dq; + + let c = self.c(&q_partial); + + result[index] = (c - c_original) / dq + } + } + + result + } + + fn compute( + &mut self, + q: &DVector, + q_prev: &DVector, + dt: Scalar, + ) -> (Scalar, Scalar, RowDVector, RowDVector) { + let c = self.c(q); + let c_prev = self.c(q_prev); + + let c_dot = (c - c_prev) / dt; + + let jacobian = self.partial_derivative(&q); + + let jacobian_dot = (jacobian.clone() - self.jacobian_prev()) / dt; + self.set_jacobian(jacobian.clone()); + + (c, c_dot, jacobian, jacobian_dot) + } +} + +impl ParticleSystem { + pub fn enforce_constraints(&mut self, dt: Scalar) { + if self.constraints.len() == 0 { + return + } + + let size = self.particles.len() * N; + let mut q = DVector::zeros(size); + let mut q_prev = DVector::zeros(size); + + for (p, particle) in self.particles.iter().enumerate() { + for i in 0..N { + q[p * N + i] = particle.position[i]; + q_prev[p * N + i] = particle.position[i] - particle.velocity[i] * dt; + } + } + + let q_dot = (q.clone() - q_prev.clone()) / dt; + + let mut c = DVector::zeros(self.constraints.len()); + let mut c_dot = DVector::zeros(self.constraints.len()); + let mut jacobian = DMatrix::::zeros(self.constraints.len(), size); + let mut jacobian_dot = DMatrix::::zeros(self.constraints.len(), size); + + for (constraint_id, constraint) in self.constraints.iter_mut().enumerate() { + let (value, derivative, j, jdot) = constraint.compute(&q, &q_prev, dt); + jacobian.set_row(constraint_id, &j); + jacobian_dot.set_row(constraint_id, &jdot); + c[constraint_id] = value; + c_dot[constraint_id] = derivative; + } + + // println!("C {}\nC' {}", c, c_dot); + // println!("J = {}", jacobian.clone()); + // println!("J' = {}", jacobian_dot.clone()); + + let mut w = DMatrix::zeros(size, size); + for (i, particle) in self.particles.iter().enumerate() { + for j in 0..N { + *w.index_mut((i * N + j, i * N + j)) = 1.0 / particle.mass; + } + } + + let mut forces = DVector::zeros(size); + for (p, particle) in self.particles.iter().enumerate() { + for i in 0..N { + forces[p * N + i] = particle.force[i]; + } + } + + let lhs = jacobian.clone() * w.clone() * jacobian.transpose(); + // println!("LHS {}", lhs); + + let rhs = -(jacobian_dot * q_dot + + jacobian.clone() * w * forces + + SPRING_CONSTANT * c + + DAMPING_CONSTANT * c_dot); + + // println!("RHS {}", rhs); + + match lhs.lu().solve(&rhs) { + Some(lambda) => { + // println!("Lambda = {}", lambda); + let constraint_force = jacobian.transpose() * lambda; + + for i in 0..self.particles.len() { + let mut force = Vector::zeros(); + for j in 0..N { + force[j] = constraint_force[i * N + j]; + } + + self.particles[i].apply_force(force); + } + } + None => println!("Lambda not found"), + } + } +} diff --git a/physics/src/constraint/slider.rs b/physics/src/constraint/slider.rs new file mode 100644 index 0000000..fa7e840 --- /dev/null +++ b/physics/src/constraint/slider.rs @@ -0,0 +1,47 @@ +use nalgebra::{DVector, RowDVector}; + +use crate::{ + algebra::{subspace::Line, Point, Scalar, Vector, N}, + particle_system::ParticleSystem, +}; + +use super::Constraint; + +pub struct SliderConstraint { + pub particle_id: usize, + + pub line: Line, + + jacobian: RowDVector, +} + +impl ParticleSystem { + pub fn add_slider_constraint(&mut self, particle_id: usize, direction: Vector) { + let point = self.particles[particle_id].position; + self.constraints.push(Box::new(SliderConstraint { + particle_id, + line: Line::new(point, [direction]), + jacobian: RowDVector::zeros(self.particles.len() * N), + })); + } +} + +impl Constraint for SliderConstraint { + fn get_particles(&self) -> Vec { + vec![self.particle_id] + } + + fn c(&self, q: &DVector) -> Scalar { + let position = q.fixed_rows::(self.particle_id * N); + let point = Point::from_coordinates(position.into()); + self.line.distance(point) + } + + fn set_jacobian(&mut self, jacobian: RowDVector) { + self.jacobian = jacobian + } + + fn jacobian_prev(&self) -> RowDVector { + self.jacobian.clone() + } +} diff --git a/physics/src/force/drag.rs b/physics/src/force/drag.rs new file mode 100644 index 0000000..c54dc40 --- /dev/null +++ b/physics/src/force/drag.rs @@ -0,0 +1,15 @@ +use crate::algebra::Scalar; + +use super::Force; + +pub struct Drag { + pub coefficient: Scalar, +} + +impl Force for Drag { + fn apply(&self, particles: &mut Vec) { + for particle in particles { + particle.apply_force(-self.coefficient * particle.velocity); + } + } +} diff --git a/physics/src/force/gravity.rs b/physics/src/force/gravity.rs new file mode 100644 index 0000000..dd600a4 --- /dev/null +++ b/physics/src/force/gravity.rs @@ -0,0 +1,15 @@ +use crate::algebra::Vector; + +use super::Force; + +pub struct Gravity { + pub vector: Vector, +} + +impl Force for Gravity { + fn apply(&self, particles: &mut Vec) { + for particle in particles { + particle.apply_force(self.vector * particle.mass); + } + } +} diff --git a/physics/src/force/mod.rs b/physics/src/force/mod.rs new file mode 100644 index 0000000..ce10f9f --- /dev/null +++ b/physics/src/force/mod.rs @@ -0,0 +1,24 @@ +use crate::particle_system::{Particle, ParticleSystem}; + +pub mod gravity; +pub mod drag; +pub mod spring; + +pub trait Force { + fn apply(&self, particles: &mut Vec); +} + +impl ParticleSystem { + fn reset_forces(&mut self) { + for particle in &mut self.particles { + particle.reset_force(); + } + } + pub fn apply_forces(&mut self) { + self.reset_forces(); + + for force in &self.forces { + force.apply(&mut self.particles) + } + } +} diff --git a/physics/src/force/spring.rs b/physics/src/force/spring.rs new file mode 100644 index 0000000..6c2cf67 --- /dev/null +++ b/physics/src/force/spring.rs @@ -0,0 +1,28 @@ +use crate::algebra::Scalar; + +use super::Force; + +pub struct Spring { + pub particle_ids: [usize; 2], + pub rest_length: Scalar, + pub spring_constant: Scalar, + pub damping_constant: Scalar, +} + +impl Force for Spring { + fn apply(&self, particles: &mut Vec) { + let a = &particles[self.particle_ids[0]]; + let b = &particles[self.particle_ids[1]]; + let i = a.position - b.position; + let i_dot = a.velocity - b.velocity; + let i_norm = i.norm(); + + let force = -(self.spring_constant * (i_norm - self.rest_length) + + (self.damping_constant * i.dot(&i_dot) / i_norm)) + * i + / i_norm; + + particles[self.particle_ids[0]].apply_force(force); + particles[self.particle_ids[1]].apply_force(-force); + } +} diff --git a/physics/src/lib.rs b/physics/src/lib.rs new file mode 100644 index 0000000..8c24a10 --- /dev/null +++ b/physics/src/lib.rs @@ -0,0 +1,9 @@ +pub use nalgebra; + +pub mod algebra; +pub mod constraint; +pub mod force; +pub mod particle_system; +mod ppm; +pub mod renderer; +pub mod solver; diff --git a/physics/src/particle_system.rs b/physics/src/particle_system.rs new file mode 100644 index 0000000..30255ec --- /dev/null +++ b/physics/src/particle_system.rs @@ -0,0 +1,53 @@ +use crate::{ + algebra::{Point, Scalar, Vector}, + constraint::Constraint, + force::Force, +}; + +#[derive(Debug)] +pub struct Particle { + pub mass: Scalar, + pub position: Point, + pub velocity: Vector, + + /// Force accumulator + pub force: Vector, +} + +impl Particle { + pub fn new(position: Point, mass: Scalar) -> Self { + Self { + mass, + position, + velocity: Vector::zeros(), + force: Vector::zeros(), + } + } + + pub fn apply_force(&mut self, force: Vector) { + self.force += force; + } + pub fn reset_force(&mut self) { + self.force = Vector::zeros() + } +} + +// #[derive(Debug)] +pub struct ParticleSystem { + pub particles: Vec, + pub constraints: Vec>, + pub forces: Vec>, + + /// Simulation clock + pub t: Scalar, +} + +impl ParticleSystem { + pub fn get_kinetic_energy(&self) -> Scalar { + self.particles.iter().fold(0.0, |acc, p| { + let velocity = p.velocity.norm(); + let energy = p.mass * velocity * velocity / 2.0; + acc + energy + }) + } +} diff --git a/physics/src/ppm.rs b/physics/src/ppm.rs new file mode 100644 index 0000000..be26127 --- /dev/null +++ b/physics/src/ppm.rs @@ -0,0 +1,67 @@ +use std::{fs::File, io::Write, path::PathBuf}; + +use nalgebra::Point2; +use crate::algebra::{Scalar, N}; +use crate::renderer::Camera; +use crate::particle_system::Particle; + +pub struct PPM { + pub prefix: PathBuf, + + buffer: [[bool; WIDTH]; HEIGHT], + // pixels_per_unit: usize, +} + +impl PPM { + pub fn new(prefix: PathBuf) -> Self { + Self { + prefix, + buffer: [[false; WIDTH]; HEIGHT], + } + } + + pub fn clear_buffer(&mut self) { + self.buffer = [[false; WIDTH]; HEIGHT] + } + + pub fn draw_circle(&mut self, center: Point2, radius: Scalar) { + let screen_center = Point2::new((WIDTH / 2) as Scalar, (HEIGHT / 2) as Scalar); + + for pixel_row in 0..HEIGHT { + for pixel_col in 0..WIDTH { + let point = Point2::new(pixel_col as Scalar, pixel_row as Scalar); + + if (point - center - screen_center.coords).norm() <= radius { + self.buffer[HEIGHT - pixel_row - 1][pixel_col] = true + } + } + } + } + + pub fn render_particles(&mut self, particles: &Vec, camera: &Camera) { + for p in particles { + let point = camera.world_to_screen_space(p.position); + self.draw_circle(point, p.mass.powf(1.0 / N as Scalar)); + } + } + + pub fn save_frame(&self, time: Scalar) { + let file_name = format!("frame-{:08.3}", time); + let path = self.prefix.join(file_name); + let mut file = File::create(path).unwrap(); + + let mut s = format!("P3\n{} {}\n255\n", WIDTH, HEIGHT); + let white = "255 255 255 "; + let black = "0 0 0 "; + + for row in self.buffer { + for pixel in row { + let color = if pixel { black } else { white }; + s += color + } + s += "\n"; + } + + file.write(s.as_bytes()).unwrap(); + } +} diff --git a/physics/src/renderer/mod.rs b/physics/src/renderer/mod.rs new file mode 100644 index 0000000..63894b3 --- /dev/null +++ b/physics/src/renderer/mod.rs @@ -0,0 +1,68 @@ +use nalgebra::{Point2, SMatrix}; + +use crate::{ + algebra::subspace::Plane, + algebra::{Point, Scalar, Vector, N}, +}; + +pub struct Camera { + plane: Plane, + up: Vector, + origin: Point, + world_to_screen_space: SMatrix, + screen_space_to_world: SMatrix, +} + +impl Camera { + pub fn new(origin: Point, up: Vector, right: Vector) -> Self { + assert!( + up.dot(&right) == 0.0, + "Up and right vectors must be orthogonal" + ); + let plane = Plane::new(origin, [up, right]); + + let screen_space_to_world = SMatrix::::from_columns(&[right, up]); + let world_to_screen_space = screen_space_to_world.pseudo_inverse(0.001).unwrap(); + + Self { + plane, + up, + origin, + world_to_screen_space, + screen_space_to_world, + } + } + + pub fn world_to_screen_space(&self, point: Point) -> Point2 { + let projected = self.plane.project_point(point); + let in_screen_space = self.world_to_screen_space * (projected - self.origin); + in_screen_space.into() + } + pub fn screen_space_to_world(&self, point: Point2) -> Point { + (self.screen_space_to_world * point) + self.origin.coords + } +} + +#[cfg(test)] +mod tests { + use crate::algebra::{Point, Vector}; + + use super::Camera; + + #[test] + fn test_projection() { + let camera = Camera::new( + Point::new(1.0, 0.0), + Vector::new(1.0, 2.0), + Vector::new(2.0, -1.0), + ); + + let point = Point::new(3.0, 1.0); + + let diff = camera.world_to_screen_space * point - Point::new(1.0, 1.0); + assert!( + diff.norm() < 0.001, + "Camera translated point into screen_space incorrectly" + ); + } +} diff --git a/physics/src/solver/midpoint.rs b/physics/src/solver/midpoint.rs new file mode 100644 index 0000000..2d71758 --- /dev/null +++ b/physics/src/solver/midpoint.rs @@ -0,0 +1,18 @@ +use crate::{algebra::Scalar, particle_system::ParticleSystem}; +use super::{PhaseSpace, Solver}; + +impl Solver for ParticleSystem { + fn step(&mut self, dt: Scalar) { + let mut state = self.collect_phase_space(); + + // Shift to the midpoint + self.scatter_phase_space(&PhaseSpace { + 0: state.0.clone() + self.compute_derivative().0 * dt / 2.0, + }); + + state.0 += self.compute_derivative().0 * dt; + self.scatter_phase_space(&state); + + self.t += dt; + } +} diff --git a/physics/src/solver/mod.rs b/physics/src/solver/mod.rs new file mode 100644 index 0000000..726dcae --- /dev/null +++ b/physics/src/solver/mod.rs @@ -0,0 +1,188 @@ +use crate::particle_system::ParticleSystem; +use crate::algebra::{Point, Scalar, Vector, N}; +use nalgebra::{Const, DVector, Dyn, Matrix, ViewStorage}; + +mod midpoint; + +/// A vector of concatenated position and velocity components of each particle +#[derive(Debug, Clone)] +pub struct PhaseSpace(DVector); +type ParticleView<'a> = Matrix< + Scalar, + Const<{ PhaseSpace::PARTICLE_DIM }>, + Const<1>, + ViewStorage<'a, Scalar, Const<{ PhaseSpace::PARTICLE_DIM }>, Const<1>, Const<1>, Dyn>, +>; + +impl PhaseSpace { + /// Each particle spans 2N elements in a vector + /// first N for position, then N more for velocity + const PARTICLE_DIM: usize = N * 2; + + pub fn new(particle_count: usize) -> Self { + let dimension = particle_count * PhaseSpace::PARTICLE_DIM; + Self(DVector::::zeros(dimension)) + } + + pub fn particle_view(&self, i: usize) -> ParticleView { + self.0 + .fixed_rows::<{ PhaseSpace::PARTICLE_DIM }>(i * PhaseSpace::PARTICLE_DIM) + } + + pub fn set_particle(&mut self, i: usize, position: Point, velocity: Vector) { + let mut view = self + .0 + .fixed_rows_mut::<{ PhaseSpace::PARTICLE_DIM }>(i * PhaseSpace::PARTICLE_DIM); + for i in 0..N { + view[i] = position[i]; + view[i + N] = velocity[i]; + } + } +} + +impl ParticleSystem { + fn collect_phase_space(&self) -> PhaseSpace { + let mut phase_space = PhaseSpace::new(self.particles.len()); + for (particle_index, particle) in self.particles.iter().enumerate() { + phase_space.set_particle(particle_index, particle.position, particle.velocity); + } + phase_space + } + + fn compute_derivative(&self) -> PhaseSpace { + let mut phase_space = PhaseSpace::new(self.particles.len()); + for (particle_index, particle) in self.particles.iter().enumerate() { + phase_space.set_particle( + particle_index, + particle.velocity.into(), + particle.force / particle.mass, + ); + } + phase_space + } + + fn scatter_phase_space(&mut self, phase_space: &PhaseSpace) { + for (particle_index, particle) in &mut self.particles.iter_mut().enumerate() { + let view = phase_space.particle_view(particle_index); + for i in 0..N { + particle.position[i] = view[i]; + particle.velocity[i] = view[i + N]; + } + } + } +} + +pub trait Solver { + fn step(&mut self, dt: Scalar); +} + +#[cfg(test)] +mod tests { + use super::{ParticleSystem, PhaseSpace, Point, Scalar, Solver, Vector}; + use crate::particle_system::Particle; + + #[test] + fn test_collect_phase_space() { + let system = ParticleSystem { + particles: vec![Particle::new(Point::new(2.0, 3.0), 1.0)], + constraints: vec![], + forces: vec![], + t: 0.0, + }; + let phase_space = system.collect_phase_space(); + + assert!( + !phase_space.0.is_empty(), + "Phase space has to contain non-zero values" + ); + } + + #[test] + fn test_scatter_phase_space() { + let mut phase_space = PhaseSpace::new(2); + phase_space.set_particle(1, Point::new(5.0, 7.0), Vector::x()); + + let mut system = ParticleSystem { + particles: vec![ + Particle::new(Point::new(0.0, 0.0), 1.0), + Particle::new(Point::new(0.0, 0.0), 1.0), + ], + constraints: vec![], + forces: vec![], + t: 0.0, + }; + + system.scatter_phase_space(&phase_space); + + assert!( + !system.particles[1].velocity.is_empty(), + "Velocity has to be set" + ); + assert!( + !system.particles[1].position.is_empty(), + "Position has to be set" + ); + } + + fn simulate_falling_ball( + fall_time: Scalar, + dt: Scalar, + mass: Scalar, + ) -> (Vector, Vector, ParticleSystem) { + let gravity = -9.8 * Vector::y(); + + let mut system = ParticleSystem { + particles: vec![Particle::new(Point::origin(), mass)], + constraints: vec![], + forces: vec![], + t: 0.0, + }; + + let iterations = (fall_time / dt) as usize; + + for _ in 0..iterations { + for particle in &mut system.particles { + particle.reset_force(); + particle.apply_force(gravity * particle.mass); + } + system.step(dt); + } + + let expected_velocity = gravity * fall_time; // vt + let expected_position = gravity * fall_time * fall_time / 2.0; // at^2 / 2 + + ( + system.particles[0].position.coords - expected_position, + system.particles[0].velocity - expected_velocity, + system, + ) + } + + #[test] + fn ball_should_fall() { + let (position_error, velocity_error, _) = simulate_falling_ball(10.0, 0.01, 3.0); + assert!( + position_error.norm() < 0.01, + "Position error is too high: {}", + position_error, + ); + assert!( + velocity_error.norm() < 0.01, + "Velocity error is too high: {}", + velocity_error, + ); + } + + #[test] + fn freefall_different_masses() { + let (_, _, system1) = simulate_falling_ball(10.0, 0.01, 2.0); + let (_, _, system2) = simulate_falling_ball(10.0, 0.01, 10.0); + + let diff = system1.particles[0].position - system2.particles[0].position; + assert!( + diff.norm() < 0.01, + "Points with different masses should fall with the same speed, diff: {}", + diff + ); + } +} diff --git a/playground/.gitignore b/playground/.gitignore new file mode 100644 index 0000000..c868ffa --- /dev/null +++ b/playground/.gitignore @@ -0,0 +1,2 @@ +/target +out diff --git a/playground/Cargo.lock b/playground/Cargo.lock new file mode 100644 index 0000000..51c385d --- /dev/null +++ b/playground/Cargo.lock @@ -0,0 +1,783 @@ +# This file is automatically @generated by Cargo. +# It is not intended for manual editing. +version = 3 + +[[package]] +name = "ahash" +version = "0.3.8" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e8fd72866655d1904d6b0997d0b07ba561047d070fbe29de039031c641b61217" + +[[package]] +name = "aho-corasick" +version = "1.1.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8e60d3430d3a69478ad0993f19238d2df97c507009a52b3c10addcd7f6bcb916" +dependencies = [ + "memchr", +] + +[[package]] +name = "approx" +version = "0.5.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "cab112f0a86d568ea0e627cc1d6be74a1e9cd55214684db5561995f6dad897c6" +dependencies = [ + "num-traits", +] + +[[package]] +name = "arrayvec" +version = "0.5.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "23b62fc65de8e4e7f52534fb52b0f3ed04746ae267519eef2a83941e8085068b" + +[[package]] +name = "autocfg" +version = "1.4.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ace50bade8e6234aa140d9a2f552bbee1db4d353f69b8217bc503490fc1a9f26" + +[[package]] +name = "bindgen" +version = "0.69.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "271383c67ccabffb7381723dea0672a673f292304fcb45c01cc648c7a8d58088" +dependencies = [ + "bitflags", + "cexpr", + "clang-sys", + "itertools", + "lazy_static", + "lazycell", + "log", + "prettyplease", + "proc-macro2", + "quote", + "regex", + "rustc-hash", + "shlex", + "syn 2.0.90", + "which", +] + +[[package]] +name = "bitflags" +version = "2.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b048fb63fd8b5923fc5aa7b340d8e156aec7ec02f0c78fa8a6ddc2613f6f71de" + +[[package]] +name = "bytemuck" +version = "1.20.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8b37c88a63ffd85d15b406896cc343916d7cf57838a847b3a6f2ca5d39a5695a" + +[[package]] +name = "cc" +version = "1.2.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9157bbaa6b165880c27a4293a474c91cdcf265cc68cc829bf10be0964a391caf" +dependencies = [ + "shlex", +] + +[[package]] +name = "cexpr" +version = "0.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6fac387a98bb7c37292057cffc56d62ecb629900026402633ae9160df93a8766" +dependencies = [ + "nom", +] + +[[package]] +name = "cfg-if" +version = "0.1.10" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4785bdd1c96b2a846b2bd7cc02e86b6b3dbf14e7e53446c4f54c92a361040822" + +[[package]] +name = "cfg-if" +version = "1.0.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "baf1de4339761588bc0619e3cbc0120ee582ebb74b53b4efbf79117bd2da40fd" + +[[package]] +name = "clang-sys" +version = "1.8.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0b023947811758c97c59bf9d1c188fd619ad4718dcaa767947df1cadb14f39f4" +dependencies = [ + "glob", + "libc", + "libloading", +] + +[[package]] +name = "cmake" +version = "0.1.52" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c682c223677e0e5b6b7f63a64b9351844c3f1b1678a68b7ee617e30fb082620e" +dependencies = [ + "cc", +] + +[[package]] +name = "crossbeam-queue" +version = "0.2.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "774ba60a54c213d409d5353bda12d49cd68d14e45036a285234c8d6f91f92570" +dependencies = [ + "cfg-if 0.1.10", + "crossbeam-utils", + "maybe-uninit", +] + +[[package]] +name = "crossbeam-utils" +version = "0.7.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c3c7c73a2d1e9fc0886a08b93e98eb643461230d5f1925e4036204d5f2e261a8" +dependencies = [ + "autocfg", + "cfg-if 0.1.10", + "lazy_static", +] + +[[package]] +name = "either" +version = "1.13.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "60b1af1c220855b6ceac025d3f6ecdd2b7c4894bfe9cd9bda4fbb4bc7c0d4cf0" + +[[package]] +name = "errno" +version = "0.3.10" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "33d852cb9b869c2a9b3df2f71a3074817f01e1844f839a144f5fcef059a4eb5d" +dependencies = [ + "libc", + "windows-sys 0.59.0", +] + +[[package]] +name = "fs_extra" +version = "1.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "42703706b716c37f96a77aea830392ad231f44c9e9a67872fa5548707e11b11c" + +[[package]] +name = "glob" +version = "0.3.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d2fabcfbdc87f4758337ca535fb41a6d701b65693ce38287d856d1674551ec9b" + +[[package]] +name = "hashbrown" +version = "0.7.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "96282e96bfcd3da0d3aa9938bedf1e50df3269b6db08b4876d2da0bb1a0841cf" +dependencies = [ + "ahash", + "autocfg", +] + +[[package]] +name = "hibitset" +version = "0.6.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f3ede5cfa60c958e60330d65163adbc4211e15a2653ad80eb0cce878de120121" + +[[package]] +name = "home" +version = "0.5.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e3d1354bf6b7235cb4a0576c2619fd4ed18183f689b12b006a0ee7329eeff9a5" +dependencies = [ + "windows-sys 0.52.0", +] + +[[package]] +name = "itertools" +version = "0.12.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ba291022dbbd398a455acf126c1e341954079855bc60dfdda641363bd6922569" +dependencies = [ + "either", +] + +[[package]] +name = "lazy_static" +version = "1.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bbd2bcb4c963f2ddae06a2efc7e9f3591312473c50c6685e1f298068316e66fe" + +[[package]] +name = "lazycell" +version = "1.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "830d08ce1d1d941e6b30645f1a0eb5643013d835ce3779a5fc208261dbe10f55" + +[[package]] +name = "libc" +version = "0.2.168" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5aaeb2981e0606ca11d79718f8bb01164f1d6ed75080182d3abf017e6d244b6d" + +[[package]] +name = "libloading" +version = "0.8.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "fc2f4eb4bc735547cfed7c0a4922cbd04a4655978c09b54f1f7b228750664c34" +dependencies = [ + "cfg-if 1.0.0", + "windows-targets", +] + +[[package]] +name = "linux-raw-sys" +version = "0.4.14" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "78b3ae25bc7c8c38cec158d1f2757ee79e9b3740fbc7ccf0e59e4b08d793fa89" + +[[package]] +name = "lock_api" +version = "0.4.12" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "07af8b9cdd281b7915f413fa73f29ebd5d55d0d3f0155584dade1ff18cea1b17" +dependencies = [ + "autocfg", + "scopeguard", +] + +[[package]] +name = "log" +version = "0.4.22" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a7a70ba024b9dc04c27ea2f0c0548feb474ec5c54bba33a7f72f873a39d07b24" + +[[package]] +name = "matrixmultiply" +version = "0.3.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9380b911e3e96d10c1f415da0876389aaf1b56759054eeb0de7df940c456ba1a" +dependencies = [ + "autocfg", + "rawpointer", +] + +[[package]] +name = "maybe-uninit" +version = "2.0.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "60302e4db3a61da70c0cb7991976248362f30319e88850c487b9b95bbf059e00" + +[[package]] +name = "memchr" +version = "2.7.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "78ca9ab1a0babb1e7d5695e3530886289c18cf2f87ec19a575a0abdce112e3a3" + +[[package]] +name = "minimal-lexical" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "68354c5c6bd36d73ff3feceb05efa59b6acb7626617f4962be322a825e61f79a" + +[[package]] +name = "mopa" +version = "0.2.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a785740271256c230f57462d3b83e52f998433a7062fc18f96d5999474a9f915" + +[[package]] +name = "nalgebra" +version = "0.33.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "26aecdf64b707efd1310e3544d709c5c0ac61c13756046aaaba41be5c4f66a3b" +dependencies = [ + "approx", + "matrixmultiply", + "nalgebra-macros", + "num-complex", + "num-rational", + "num-traits", + "simba", + "typenum", +] + +[[package]] +name = "nalgebra-macros" +version = "0.2.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "254a5372af8fc138e36684761d3c0cdb758a4410e938babcff1c860ce14ddbfc" +dependencies = [ + "proc-macro2", + "quote", + "syn 2.0.90", +] + +[[package]] +name = "nom" +version = "7.1.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d273983c5a657a70a3e8f2a01329822f3b8c8172b73826411a55751e404a0a4a" +dependencies = [ + "memchr", + "minimal-lexical", +] + +[[package]] +name = "num-bigint" +version = "0.4.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a5e44f723f1133c9deac646763579fdb3ac745e418f2a7af9cd0c431da1f20b9" +dependencies = [ + "num-integer", + "num-traits", +] + +[[package]] +name = "num-complex" +version = "0.4.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "73f88a1307638156682bada9d7604135552957b7818057dcef22705b4d509495" +dependencies = [ + "num-traits", +] + +[[package]] +name = "num-integer" +version = "0.1.46" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7969661fd2958a5cb096e56c8e1ad0444ac2bbcd0061bd28660485a44879858f" +dependencies = [ + "num-traits", +] + +[[package]] +name = "num-rational" +version = "0.4.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f83d14da390562dca69fc84082e73e548e1ad308d24accdedd2720017cb37824" +dependencies = [ + "num-bigint", + "num-integer", + "num-traits", +] + +[[package]] +name = "num-traits" +version = "0.2.19" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "071dfc062690e90b734c0b2273ce72ad0ffa95f0c74596bc250dcfd960262841" +dependencies = [ + "autocfg", +] + +[[package]] +name = "once_cell" +version = "1.20.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1261fe7e33c73b354eab43b1273a57c8f967d0391e80353e51f764ac02cf6775" + +[[package]] +name = "parking_lot" +version = "0.12.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f1bf18183cf54e8d6059647fc3063646a1801cf30896933ec2311622cc4b9a27" +dependencies = [ + "lock_api", + "parking_lot_core", +] + +[[package]] +name = "parking_lot_core" +version = "0.9.10" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1e401f977ab385c9e4e3ab30627d6f26d00e2c73eef317493c4ec6d468726cf8" +dependencies = [ + "cfg-if 1.0.0", + "libc", + "redox_syscall", + "smallvec", + "windows-targets", +] + +[[package]] +name = "paste" +version = "1.0.15" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "57c0d7b74b563b49d38dae00a0c37d4d6de9b432382b2892f0574ddcae73fd0a" + +[[package]] +name = "physics" +version = "0.1.0" +dependencies = [ + "nalgebra", +] + +[[package]] +name = "playground" +version = "0.1.0" +dependencies = [ + "physics", + "raylib", +] + +[[package]] +name = "prettyplease" +version = "0.2.25" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "64d1ec885c64d0457d564db4ec299b2dae3f9c02808b8ad9c3a089c591b18033" +dependencies = [ + "proc-macro2", + "syn 2.0.90", +] + +[[package]] +name = "proc-macro2" +version = "1.0.92" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "37d3544b3f2748c54e147655edb5025752e2303145b5aefb3c3ea2c78b973bb0" +dependencies = [ + "unicode-ident", +] + +[[package]] +name = "quote" +version = "1.0.37" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b5b9d34b8991d19d98081b46eacdd8eb58c6f2b201139f7c5f643cc155a633af" +dependencies = [ + "proc-macro2", +] + +[[package]] +name = "rawpointer" +version = "0.2.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "60a357793950651c4ed0f3f52338f53b2f809f32d83a07f72909fa13e4c6c1e3" + +[[package]] +name = "raylib" +version = "5.0.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e2a7a6734329d7b872a418fe4cb08ca282eb66a6f4a3430bd4ee4e6a8cac6632" +dependencies = [ + "cfg-if 1.0.0", + "lazy_static", + "libc", + "parking_lot", + "raylib-sys", + "specs", + "specs-derive", +] + +[[package]] +name = "raylib-sys" +version = "5.0.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "db5c6001cfaeec17210713227d11f3b1ba4b723bb12cff47d1b93c4060e10ad0" +dependencies = [ + "bindgen", + "cc", + "cmake", + "fs_extra", +] + +[[package]] +name = "redox_syscall" +version = "0.5.8" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "03a862b389f93e68874fbf580b9de08dd02facb9a788ebadaf4a3fd33cf58834" +dependencies = [ + "bitflags", +] + +[[package]] +name = "regex" +version = "1.11.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b544ef1b4eac5dc2db33ea63606ae9ffcfac26c1416a2806ae0bf5f56b201191" +dependencies = [ + "aho-corasick", + "memchr", + "regex-automata", + "regex-syntax", +] + +[[package]] +name = "regex-automata" +version = "0.4.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "809e8dc61f6de73b46c85f4c96486310fe304c434cfa43669d7b40f711150908" +dependencies = [ + "aho-corasick", + "memchr", + "regex-syntax", +] + +[[package]] +name = "regex-syntax" +version = "0.8.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2b15c43186be67a4fd63bee50d0303afffcef381492ebe2c5d87f324e1b8815c" + +[[package]] +name = "rustc-hash" +version = "1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "08d43f7aa6b08d49f382cde6a7982047c3426db949b1424bc4b7ec9ae12c6ce2" + +[[package]] +name = "rustix" +version = "0.38.42" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f93dc38ecbab2eb790ff964bb77fa94faf256fd3e73285fd7ba0903b76bedb85" +dependencies = [ + "bitflags", + "errno", + "libc", + "linux-raw-sys", + "windows-sys 0.59.0", +] + +[[package]] +name = "safe_arch" +version = "0.7.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c3460605018fdc9612bce72735cba0d27efbcd9904780d44c7e3a9948f96148a" +dependencies = [ + "bytemuck", +] + +[[package]] +name = "scopeguard" +version = "1.2.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "94143f37725109f92c262ed2cf5e59bce7498c01bcc1502d7b9afe439a4e9f49" + +[[package]] +name = "shlex" +version = "1.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0fda2ff0d084019ba4d7c6f371c95d8fd75ce3524c3cb8fb653a3023f6323e64" + +[[package]] +name = "shred" +version = "0.10.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c5f08237e667ac94ad20f8878b5943d91a93ccb231428446c57c21c57779016d" +dependencies = [ + "arrayvec", + "hashbrown", + "mopa", + "smallvec", + "tynm", +] + +[[package]] +name = "shrev" +version = "1.1.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a5ea33232fdcf1bf691ca33450e5a94dde13e1a8cbb8caabc5e4f9d761e10b1a" + +[[package]] +name = "simba" +version = "0.9.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b3a386a501cd104797982c15ae17aafe8b9261315b5d07e3ec803f2ea26be0fa" +dependencies = [ + "approx", + "num-complex", + "num-traits", + "paste", + "wide", +] + +[[package]] +name = "smallvec" +version = "1.13.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3c5e1a9a646d36c3599cd173a41282daf47c44583ad367b8e6837255952e5c67" + +[[package]] +name = "specs" +version = "0.16.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "fff28a29366aff703d5da8a7e2c8875dc8453ac1118f842cbc0fa70c7db51240" +dependencies = [ + "crossbeam-queue", + "hashbrown", + "hibitset", + "log", + "shred", + "shrev", + "tuple_utils", +] + +[[package]] +name = "specs-derive" +version = "0.4.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3e23e09360f3d2190fec4222cd9e19d3158d5da948c0d1ea362df617dd103511" +dependencies = [ + "proc-macro2", + "quote", + "syn 1.0.109", +] + +[[package]] +name = "syn" +version = "1.0.109" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "72b64191b275b66ffe2469e8af2c1cfe3bafa67b529ead792a6d0160888b4237" +dependencies = [ + "proc-macro2", + "quote", + "unicode-ident", +] + +[[package]] +name = "syn" +version = "2.0.90" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "919d3b74a5dd0ccd15aeb8f93e7006bd9e14c295087c9896a110f490752bcf31" +dependencies = [ + "proc-macro2", + "quote", + "unicode-ident", +] + +[[package]] +name = "tuple_utils" +version = "0.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "44834418e2c5b16f47bedf35c28e148db099187dd5feee6367fb2525863af4f1" + +[[package]] +name = "tynm" +version = "0.1.10" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "bd30d05e69d1478e13fe3e7a853409cfec82cebc2cf9b8d613b3c6b0081781ed" +dependencies = [ + "nom", +] + +[[package]] +name = "typenum" +version = "1.17.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "42ff0bf0c66b8238c6f3b578df37d0b7848e55df8577b3f74f92a69acceeb825" + +[[package]] +name = "unicode-ident" +version = "1.0.14" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "adb9e6ca4f869e1180728b7950e35922a7fc6397f7b641499e8f3ef06e50dc83" + +[[package]] +name = "which" +version = "4.4.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "87ba24419a2078cd2b0f2ede2691b6c66d8e47836da3b6db8265ebad47afbfc7" +dependencies = [ + "either", + "home", + "once_cell", + "rustix", +] + +[[package]] +name = "wide" +version = "0.7.30" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "58e6db2670d2be78525979e9a5f9c69d296fd7d670549fe9ebf70f8708cb5019" +dependencies = [ + "bytemuck", + "safe_arch", +] + +[[package]] +name = "windows-sys" +version = "0.52.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "282be5f36a8ce781fad8c8ae18fa3f9beff57ec1b52cb3de0789201425d9a33d" +dependencies = [ + "windows-targets", +] + +[[package]] +name = "windows-sys" +version = "0.59.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1e38bc4d79ed67fd075bcc251a1c39b32a1776bbe92e5bef1f0bf1f8c531853b" +dependencies = [ + "windows-targets", +] + +[[package]] +name = "windows-targets" +version = "0.52.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9b724f72796e036ab90c1021d4780d4d3d648aca59e491e6b98e725b84e99973" +dependencies = [ + "windows_aarch64_gnullvm", + "windows_aarch64_msvc", + "windows_i686_gnu", + "windows_i686_gnullvm", + "windows_i686_msvc", + "windows_x86_64_gnu", + "windows_x86_64_gnullvm", + "windows_x86_64_msvc", +] + +[[package]] +name = "windows_aarch64_gnullvm" +version = "0.52.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "32a4622180e7a0ec044bb555404c800bc9fd9ec262ec147edd5989ccd0c02cd3" + +[[package]] +name = "windows_aarch64_msvc" +version = "0.52.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "09ec2a7bb152e2252b53fa7803150007879548bc709c039df7627cabbd05d469" + +[[package]] +name = "windows_i686_gnu" +version = "0.52.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8e9b5ad5ab802e97eb8e295ac6720e509ee4c243f69d781394014ebfe8bbfa0b" + +[[package]] +name = "windows_i686_gnullvm" +version = "0.52.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0eee52d38c090b3caa76c563b86c3a4bd71ef1a819287c19d586d7334ae8ed66" + +[[package]] +name = "windows_i686_msvc" +version = "0.52.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "240948bc05c5e7c6dabba28bf89d89ffce3e303022809e73deaefe4f6ec56c66" + +[[package]] +name = "windows_x86_64_gnu" +version = "0.52.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "147a5c80aabfbf0c7d901cb5895d1de30ef2907eb21fbbab29ca94c5b08b1a78" + +[[package]] +name = "windows_x86_64_gnullvm" +version = "0.52.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "24d5b23dc417412679681396f2b49f3de8c1473deb516bd34410872eff51ed0d" + +[[package]] +name = "windows_x86_64_msvc" +version = "0.52.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "589f6da84c646204747d1270a2a5661ea66ed1cced2631d546fdfb155959f9ec" diff --git a/playground/Cargo.toml b/playground/Cargo.toml new file mode 100644 index 0000000..163ad90 --- /dev/null +++ b/playground/Cargo.toml @@ -0,0 +1,8 @@ +[package] +name = "playground" +version = "0.1.0" +edition = "2021" + +[dependencies] +physics = { path = "../physics" } +raylib = "5.0.2" diff --git a/playground/src/main.rs b/playground/src/main.rs new file mode 100644 index 0000000..1cb8198 --- /dev/null +++ b/playground/src/main.rs @@ -0,0 +1,185 @@ +use raylib::prelude::*; + +use physics::force::{drag::Drag, gravity::Gravity, spring::Spring}; +use physics::nalgebra::Point as PointBase; +use physics::particle_system::{Particle, ParticleSystem}; +use physics::algebra::{Point, Scalar, Vector, N}; +use physics::renderer::Camera; +use physics::solver::Solver; + +const SCALE: i32 = 5; +fn screen_space_to_raylib(p: PointBase, d: &RaylibDrawHandle) -> PointBase { + PointBase::::new( + d.get_screen_width() / 2 + p.x as i32 * SCALE, + d.get_screen_height() / 2 - p.y as i32 * SCALE, + ) +} +fn raylib_to_screen_space(p: PointBase, d: &RaylibDrawHandle) -> PointBase { + PointBase::::new( + ((p.x - d.get_screen_width() / 2) / SCALE) as Scalar, + ((p.y - d.get_screen_height() / 2) / -SCALE) as Scalar, + ) +} + +fn main() { + let dt = 0.0001; + let mut system = ParticleSystem { + particles: vec![ + Particle::new(Point::origin(), 4.0), + Particle::new(Point::new(-30.0, 0.0), 15.0), + Particle::new(Point::new(20.0, 0.0), 30.0), + Particle::new(Point::new(5.0, 20.0), 50.0), + Particle::new(Point::origin(), 25.0), + Particle::new(Point::new(50.0, 0.0), 10.0), + Particle::new(Point::new(-100.0, -100.0), 100.0), + ], + constraints: vec![], + forces: vec![ + Box::new(Gravity { + vector: Vector::y() * -9.8, + }), + Box::new(Drag { coefficient: 0.2 }), + Box::new(Spring { + particle_ids: [4, 2], + spring_constant: 0.75, + damping_constant: 0.1, + rest_length: 20.0, + }), + ], + t: 0.0, + }; + + system.add_anchor_constraint(0); + system.add_beam_constraint([0, 2]); + system.add_beam_constraint([1, 2]); + system.add_beam_constraint([1, 3]); + system.add_beam_constraint([2, 3]); + system.add_slider_constraint(5, Vector::x()); + system.add_beam_constraint([5, 4]); + + let mut selected_particle_id = None; + + let (mut rl, thread) = raylib::init() + .size(640, 480) + .title("Physics simulation") + .resizable() + .build(); + + while !rl.window_should_close() { + let mut d = rl.begin_drawing(&thread); + + d.clear_background(Color::WHITE); + d.draw_text( + format!( + "Time: {:03.3}\n\nFPS: {}\n\nKinetic energy: {:03.0}", + system.t, + d.get_fps(), + system.get_kinetic_energy() + ) + .as_str(), + 12, + 12, + 20, + Color::BLACK, + ); + + let camera = Camera::new(Point::origin(), Vector::y(), Vector::x()); + + match selected_particle_id { + Some(particle_id) => { + // let p: &Particle = &system.particles[particle_id]; + + let mouse_point = PointBase::::new(d.get_mouse_x(), d.get_mouse_y()); + let screen_space = raylib_to_screen_space(mouse_point, &d); + let world_mouse = camera.screen_space_to_world(screen_space); + + let mouse_particle_id = system.particles.len() - 1; + system.particles[mouse_particle_id].position = world_mouse; + + system.forces.push(Box::new(Spring { + particle_ids: [mouse_particle_id, particle_id], + spring_constant: 0.99, + damping_constant: 0.99, + rest_length: 1.0, + })); + } + None => (), + } + + for _ in 0..10 { + system.apply_forces(); + system.enforce_constraints(dt); + system.step(dt); + } + + match selected_particle_id { + Some(_) => { + system.forces.pop(); + } + None => (), + } + + if d.is_mouse_button_released(MouseButton::MOUSE_BUTTON_LEFT) { + selected_particle_id = None + } + + for (particle_id, particle) in system.particles.iter().enumerate() { + let position = + screen_space_to_raylib(camera.world_to_screen_space(particle.position), &d); + let radius = particle.mass.powf(1.0 / (N as f64)) as f32 * SCALE as f32; + + if d.is_mouse_button_down(MouseButton::MOUSE_BUTTON_LEFT) + && selected_particle_id.is_none() + { + let mouse_point = + PointBase::::new(d.get_mouse_x() as f32, d.get_mouse_y() as f32); + let position = PointBase::::new(position.x as f32, position.y as f32); + if (position - mouse_point).norm() < radius { + selected_particle_id = Some(particle_id) + } + } + + d.draw_circle( + position.x, + position.y, + radius, + if selected_particle_id.is_some_and(|v| v == particle_id) { + Color::RED + } else { + Color::BLACK + }, + ); + } + for c in &system.constraints { + let particle_ids = c.get_particles(); + if particle_ids.len() == 2 { + let a = screen_space_to_raylib(system.particles[particle_ids[0]].position, &d); + let b = screen_space_to_raylib(system.particles[particle_ids[1]].position, &d); + + d.draw_line(a.x, a.y, b.x, b.y, Color::GRAY); + } + } + + { + // Hard-coded spring + let particle_ids = [4, 2]; + let a = screen_space_to_raylib(system.particles[particle_ids[0]].position, &d); + let b = screen_space_to_raylib(system.particles[particle_ids[1]].position, &d); + d.draw_line(a.x, a.y, b.x, b.y, Color::GREEN); + } + { + // Hard-coded spring + match selected_particle_id { + Some(id) => { + let mouse_particle_id = system.particles.len() - 1; + let a = screen_space_to_raylib(system.particles[id].position, &d); + let b = + screen_space_to_raylib(system.particles[mouse_particle_id].position, &d); + + d.draw_line(a.x, a.y, b.x, b.y, Color::GREEN); + } + None => (), + } + } + } +} diff --git a/src/algebra/mod.rs b/src/algebra/mod.rs deleted file mode 100644 index 13dcd6a..0000000 --- a/src/algebra/mod.rs +++ /dev/null @@ -1 +0,0 @@ -pub mod subspace; diff --git a/src/algebra/subspace.rs b/src/algebra/subspace.rs deleted file mode 100644 index 3b32ba7..0000000 --- a/src/algebra/subspace.rs +++ /dev/null @@ -1,69 +0,0 @@ -use nalgebra::SMatrix; - -use crate::{particle_system::Scalar, particle_system::N, Point, Vector}; - -type ProjectionMatrix = SMatrix; - -pub struct Subspace { - point: Point, - vectors: [Vector; DIM], - projection_matrix: ProjectionMatrix, -} - -pub type Line = Subspace<1>; -pub type Plane = Subspace<2>; - -impl Subspace { - pub fn new(point: Point, mut vectors: [Vector; DIM]) -> Self { - for vector in &mut vectors { - vector.normalize_mut(); - } - - Self { - point, - vectors, - projection_matrix: { - let v = SMatrix::::from_columns(&vectors); - let transpose = v.transpose(); - let inverse = (transpose * v).try_inverse().unwrap(); - v * inverse * transpose - } - } - } - - pub fn project_point(&self, point: Point) -> Point { - self.projection_matrix * (point - self.point.coords) + self.point.coords - } - - pub fn distance(&self, point: Point) -> Scalar { - let projected = self.project_point(point); - (projected - point).norm() - } -} - -#[cfg(test)] -mod tests { - use crate::{algebra::subspace::Line, Point, Vector}; - - #[test] - fn test_projection_onto_line() { - let line = Line::new(Point::new(1.0, 0.0), [Vector::new(3.0, 1.0)]); - let point = Point::new(3.0, 4.0); - - let matrix = line.projection_matrix; - { - let diff = matrix * matrix - matrix; - assert!(diff.norm() < 0.001, "Projection matrix squared should be equal to itself: {}, {}", matrix, matrix* matrix); - } - { - let diff = matrix - matrix.transpose(); - assert!(diff.norm() < 0.001, "Projection matrix transpose should be equal to itself: {}, {}", matrix, matrix.transpose()); - } - - { - let projected = line.project_point(point); - let diff = projected - Point::new(4.0, 1.0); - assert!(diff.norm() < 0.001, "Point projected incorrectly"); - } - } -} diff --git a/src/constraint/anchor.rs b/src/constraint/anchor.rs deleted file mode 100644 index 3ef94ff..0000000 --- a/src/constraint/anchor.rs +++ /dev/null @@ -1,42 +0,0 @@ -use nalgebra::{DVector, RowDVector}; - -use crate::particle_system::{ParticleSystem, Point, Scalar, N}; - -use super::Constraint; - -pub struct AnchorConstraint { - pub particle_id: usize, - pub anchor: Point, - - jacobian: RowDVector, -} - -impl ParticleSystem { - pub fn add_anchor_constraint(&mut self, particle_id: usize) { - let anchor = self.particles[particle_id].position; - self.constraints.push(Box::new(AnchorConstraint { - particle_id, - anchor, - jacobian: RowDVector::zeros(self.particles.len() * N), - })); - } -} - -impl Constraint for AnchorConstraint { - fn get_particles(&self) -> Vec { - vec![self.particle_id] - } - - fn c(&self, q: &DVector) -> Scalar { - let position = q.fixed_rows(self.particle_id * N); - (position - self.anchor.coords).norm() - } - - fn set_jacobian(&mut self, jacobian: RowDVector) { - self.jacobian = jacobian - } - - fn jacobian_prev(&self) -> RowDVector { - self.jacobian.clone() - } -} diff --git a/src/constraint/beam.rs b/src/constraint/beam.rs deleted file mode 100644 index c953920..0000000 --- a/src/constraint/beam.rs +++ /dev/null @@ -1,46 +0,0 @@ -use nalgebra::{DVector, RowDVector}; - -use crate::particle_system::{ParticleSystem, Scalar, N}; - -use super::Constraint; - -pub struct BeamConstraint { - pub particle_ids: [usize; 2], - pub length: Scalar, - - jacobian: RowDVector, -} - -impl ParticleSystem { - pub fn add_beam_constraint(&mut self, particle_ids: [usize; 2]) { - let a = &self.particles[particle_ids[0]]; - let b = &self.particles[particle_ids[1]]; - - self.constraints.push(Box::new(BeamConstraint { - particle_ids, - length: (a.position - b.position).norm(), - jacobian: RowDVector::zeros(self.particles.len() * N), - })); - } -} - -impl Constraint for BeamConstraint { - fn get_particles(&self) -> Vec { - Vec::from(self.particle_ids) - } - - fn c(&self, q: &DVector) -> Scalar { - let a = q.fixed_rows::(self.particle_ids[0] * N); - let b = q.fixed_rows::(self.particle_ids[1] * N); - - (a - b).norm() - self.length - } - - fn set_jacobian(&mut self, jacobian: RowDVector) { - self.jacobian = jacobian - } - - fn jacobian_prev(&self) -> RowDVector { - self.jacobian.clone() - } -} diff --git a/src/constraint/mod.rs b/src/constraint/mod.rs deleted file mode 100644 index 6028b42..0000000 --- a/src/constraint/mod.rs +++ /dev/null @@ -1,139 +0,0 @@ -use nalgebra::{DMatrix, DVector, RowDVector}; - -use crate::particle_system::{ParticleSystem, Scalar, Vector, N}; -pub mod anchor; -pub mod slider; -pub mod beam; - -const SPRING_CONSTANT: Scalar = 0.75; -const DAMPING_CONSTANT: Scalar = 0.45; - -/// SIZE is always 3xN, but this operation can't be done at compile time yet: -/// "generic parameters may not be used in const operations" -pub trait Constraint { - fn get_particles(&self) -> Vec; - - fn c(&self, q: &DVector) -> Scalar; - - fn jacobian_prev(&self) -> RowDVector; - fn set_jacobian(&mut self, jacobian: RowDVector); - - fn partial_derivative(&self, q: &DVector) -> RowDVector { - let dq = 0.001; - let c_original = self.c(&q); - let mut result = RowDVector::zeros(q.len()); - // The only non-zero components of derivative vector are - // the particles that this constraint affects - for particle_id in self.get_particles() { - for i in 0..N { - let index = particle_id * N + i; - let mut q_partial = q.clone(); - q_partial[index] += dq; - - let c = self.c(&q_partial); - - result[index] = (c - c_original) / dq - } - } - - result - } - - fn compute( - &mut self, - q: &DVector, - q_prev: &DVector, - dt: Scalar, - ) -> (Scalar, Scalar, RowDVector, RowDVector) { - let c = self.c(q); - let c_prev = self.c(q_prev); - - let c_dot = (c - c_prev) / dt; - - let jacobian = self.partial_derivative(&q); - - let jacobian_dot = (jacobian.clone() - self.jacobian_prev()) / dt; - self.set_jacobian(jacobian.clone()); - - (c, c_dot, jacobian, jacobian_dot) - } -} - -impl ParticleSystem { - pub fn enforce_constraints(&mut self, dt: Scalar) { - if self.constraints.len() == 0 { - return - } - - let size = self.particles.len() * N; - let mut q = DVector::zeros(size); - let mut q_prev = DVector::zeros(size); - - for (p, particle) in self.particles.iter().enumerate() { - for i in 0..N { - q[p * N + i] = particle.position[i]; - q_prev[p * N + i] = particle.position[i] - particle.velocity[i] * dt; - } - } - - let q_dot = (q.clone() - q_prev.clone()) / dt; - - let mut c = DVector::zeros(self.constraints.len()); - let mut c_dot = DVector::zeros(self.constraints.len()); - let mut jacobian = DMatrix::::zeros(self.constraints.len(), size); - let mut jacobian_dot = DMatrix::::zeros(self.constraints.len(), size); - - for (constraint_id, constraint) in self.constraints.iter_mut().enumerate() { - let (value, derivative, j, jdot) = constraint.compute(&q, &q_prev, dt); - jacobian.set_row(constraint_id, &j); - jacobian_dot.set_row(constraint_id, &jdot); - c[constraint_id] = value; - c_dot[constraint_id] = derivative; - } - - // println!("C {}\nC' {}", c, c_dot); - // println!("J = {}", jacobian.clone()); - // println!("J' = {}", jacobian_dot.clone()); - - let mut w = DMatrix::zeros(size, size); - for (i, particle) in self.particles.iter().enumerate() { - for j in 0..N { - *w.index_mut((i * N + j, i * N + j)) = 1.0 / particle.mass; - } - } - - let mut forces = DVector::zeros(size); - for (p, particle) in self.particles.iter().enumerate() { - for i in 0..N { - forces[p * N + i] = particle.force[i]; - } - } - - let lhs = jacobian.clone() * w.clone() * jacobian.transpose(); - // println!("LHS {}", lhs); - - let rhs = -(jacobian_dot * q_dot - + jacobian.clone() * w * forces - + SPRING_CONSTANT * c - + DAMPING_CONSTANT * c_dot); - - // println!("RHS {}", rhs); - - match lhs.lu().solve(&rhs) { - Some(lambda) => { - // println!("Lambda = {}", lambda); - let constraint_force = jacobian.transpose() * lambda; - - for i in 0..self.particles.len() { - let mut force = Vector::zeros(); - for j in 0..N { - force[j] = constraint_force[i * N + j]; - } - - self.particles[i].apply_force(force); - } - } - None => println!("Lambda not found"), - } - } -} diff --git a/src/constraint/slider.rs b/src/constraint/slider.rs deleted file mode 100644 index 9f7973b..0000000 --- a/src/constraint/slider.rs +++ /dev/null @@ -1,50 +0,0 @@ -use std::usize; - -use nalgebra::{DVector, RowDVector}; - -use crate::{ - algebra::subspace::Line, - particle_system::{Scalar, N}, - ParticleSystem, Point, Vector, -}; - -use super::Constraint; - -pub struct SliderConstraint { - pub particle_id: usize, - - pub line: Line, - - jacobian: RowDVector, -} - -impl ParticleSystem { - pub fn add_slider_constraint(&mut self, particle_id: usize, direction: Vector) { - let point = self.particles[particle_id].position; - self.constraints.push(Box::new(SliderConstraint { - particle_id, - line: Line::new(point, [direction]), - jacobian: RowDVector::zeros(self.particles.len() * N), - })); - } -} - -impl Constraint for SliderConstraint { - fn get_particles(&self) -> Vec { - vec![self.particle_id] - } - - fn c(&self, q: &DVector) -> Scalar { - let position = q.fixed_rows::(self.particle_id * N); - let point = Point::from_coordinates(position.into()); - self.line.distance(point) - } - - fn set_jacobian(&mut self, jacobian: RowDVector) { - self.jacobian = jacobian - } - - fn jacobian_prev(&self) -> RowDVector { - self.jacobian.clone() - } -} diff --git a/src/force/drag.rs b/src/force/drag.rs deleted file mode 100644 index ef859a0..0000000 --- a/src/force/drag.rs +++ /dev/null @@ -1,15 +0,0 @@ -use crate::particle_system::Scalar; - -use super::Force; - -pub struct Drag { - pub coefficient: Scalar, -} - -impl Force for Drag { - fn apply(&self, particles: &mut Vec) { - for particle in particles { - particle.apply_force(-self.coefficient * particle.velocity); - } - } -} diff --git a/src/force/gravity.rs b/src/force/gravity.rs deleted file mode 100644 index b023e63..0000000 --- a/src/force/gravity.rs +++ /dev/null @@ -1,15 +0,0 @@ -use crate::particle_system::Vector; - -use super::Force; - -pub struct Gravity { - pub vector: Vector, -} - -impl Force for Gravity { - fn apply(&self, particles: &mut Vec) { - for particle in particles { - particle.apply_force(self.vector * particle.mass); - } - } -} diff --git a/src/force/mod.rs b/src/force/mod.rs deleted file mode 100644 index ce10f9f..0000000 --- a/src/force/mod.rs +++ /dev/null @@ -1,24 +0,0 @@ -use crate::particle_system::{Particle, ParticleSystem}; - -pub mod gravity; -pub mod drag; -pub mod spring; - -pub trait Force { - fn apply(&self, particles: &mut Vec); -} - -impl ParticleSystem { - fn reset_forces(&mut self) { - for particle in &mut self.particles { - particle.reset_force(); - } - } - pub fn apply_forces(&mut self) { - self.reset_forces(); - - for force in &self.forces { - force.apply(&mut self.particles) - } - } -} diff --git a/src/main.rs b/src/main.rs deleted file mode 100644 index 0b5704b..0000000 --- a/src/main.rs +++ /dev/null @@ -1,59 +0,0 @@ -use std::path::PathBuf; - -use force::{drag::Drag, gravity::Gravity}; -use particle_system::{Particle, ParticleSystem, Point, Vector}; -use ppm::PPM; -use solver::Solver; - -mod constraint; -mod particle_system; -mod ppm; -mod solver; -mod force; -mod algebra; - -fn main() { - let ppm = PPM { - width: 100, - height: 200, - prefix: PathBuf::from("./out"), - }; - - let dt = 0.01; - let mut system = ParticleSystem { - particles: vec![ - Particle::new(Point::origin(), 4.0), - Particle::new(Point::new(-30.0, 0.0), 15.0), - Particle::new(Point::new(20.0, 0.0), 30.0), - Particle::new(Point::new(5.0, 20.0), 50.0), - ], - constraints: vec![], - forces: vec![ - Box::new(Gravity { - vector: Vector::y() * -9.8, - }), - Box::new(Drag { coefficient: 0.0 }), - ], - t: 0.0, - }; - - system.add_anchor_constraint(0); - system.add_beam_constraint([0, 2]); - system.add_beam_constraint([1, 2]); - system.add_beam_constraint([1, 3]); - system.add_beam_constraint([2, 3]); - - for i in 0..150_00 { - system.apply_forces(); - - if i % 10 == 0 { - println!("Iteration #{i}"); - println!("{:#?}", system.particles); - ppm.save_frame(&system.particles, system.t); - } - - system.enforce_constraints(dt); - - system.step(dt); - } -} diff --git a/src/particle_system.rs b/src/particle_system.rs deleted file mode 100644 index a3d7a4b..0000000 --- a/src/particle_system.rs +++ /dev/null @@ -1,47 +0,0 @@ -use nalgebra::{Point as PointBase, SVector}; - -use crate::{constraint::Constraint, force::Force}; - -pub const N: usize = 2; -pub type Scalar = f64; - -pub type Vector = SVector; -pub type Point = PointBase; - -#[derive(Debug)] -pub struct Particle { - pub mass: Scalar, - pub position: Point, - pub velocity: Vector, - - /// Force accumulator - pub force: Vector, -} - -impl Particle { - pub fn new(position: Point, mass: Scalar) -> Self { - Self { - mass, - position, - velocity: Vector::zeros(), - force: Vector::zeros(), - } - } - - pub fn apply_force(&mut self, force: Vector) { - self.force += force; - } - pub fn reset_force(&mut self) { - self.force = Vector::zeros() - } -} - -// #[derive(Debug)] -pub struct ParticleSystem { - pub particles: Vec, - pub constraints: Vec>, - pub forces: Vec>, - - /// Simulation clock - pub t: Scalar, -} diff --git a/src/ppm.rs b/src/ppm.rs deleted file mode 100644 index 01cdebb..0000000 --- a/src/ppm.rs +++ /dev/null @@ -1,41 +0,0 @@ -use std::{fs::File, io::Write, path::PathBuf}; - -use crate::particle_system::{Particle, Vector, N, Scalar}; - -pub struct PPM { - pub prefix: PathBuf, - pub width: usize, - pub height: usize, - // pixels_per_unit: usize, -} - -impl PPM { - fn render_particles(&self, particles: &Vec) -> String { - let mut s = format!("P3\n{} {}\n255\n", self.width, self.height); - let white = "255 255 255 "; - let black = "0 0 0 "; - - for pixel_row in 0..self.height { - for pixel_col in 0..self.width { - let point = Vector::new(pixel_col as Scalar, (pixel_row as Scalar) * -1.0) - + Vector::new(self.width as Scalar / -2.0, self.height as Scalar / 2.0); - let color = match particles.iter().any(|p| { - (p.position - point).coords.norm() <= (p.mass).powf(1.0 / (N as f64)) - }) { - true => black, - false => white, - }; - s += color; - } - } - s - } - - pub fn save_frame(&self, particles: &Vec, time: Scalar) { - let file_name = format!("frame-{:08.3}", time); - let path = self.prefix.join(file_name); - let mut file = File::create(path).unwrap(); - let data = self.render_particles(particles); - file.write(data.as_bytes()).unwrap(); - } -} diff --git a/src/solver/midpoint.rs b/src/solver/midpoint.rs deleted file mode 100644 index 08a3e3c..0000000 --- a/src/solver/midpoint.rs +++ /dev/null @@ -1,18 +0,0 @@ -use crate::particle_system::{ParticleSystem, Scalar}; -use super::{PhaseSpace, Solver}; - -impl Solver for ParticleSystem { - fn step(&mut self, dt: Scalar) { - let mut state = self.collect_phase_space(); - - // Shift to the midpoint - self.scatter_phase_space(&PhaseSpace { - 0: state.0.clone() + self.compute_derivative().0 * dt / 2.0, - }); - - state.0 += self.compute_derivative().0 * dt; - self.scatter_phase_space(&state); - - self.t += dt; - } -} diff --git a/src/solver/mod.rs b/src/solver/mod.rs deleted file mode 100644 index 1544378..0000000 --- a/src/solver/mod.rs +++ /dev/null @@ -1,187 +0,0 @@ -use crate::particle_system::{ParticleSystem, Point, Scalar, Vector, N}; -use nalgebra::{Const, DVector, Dyn, Matrix, ViewStorage}; - -mod midpoint; - -/// A vector of concatenated position and velocity components of each particle -#[derive(Debug, Clone)] -pub struct PhaseSpace(DVector); -type ParticleView<'a> = Matrix< - Scalar, - Const<{ PhaseSpace::PARTICLE_DIM }>, - Const<1>, - ViewStorage<'a, Scalar, Const<{ PhaseSpace::PARTICLE_DIM }>, Const<1>, Const<1>, Dyn>, ->; - -impl PhaseSpace { - /// Each particle spans 2N elements in a vector - /// first N for position, then N more for velocity - const PARTICLE_DIM: usize = N * 2; - - pub fn new(particle_count: usize) -> Self { - let dimension = particle_count * PhaseSpace::PARTICLE_DIM; - Self(DVector::::zeros(dimension)) - } - - pub fn particle_view(&self, i: usize) -> ParticleView { - self.0 - .fixed_rows::<{ PhaseSpace::PARTICLE_DIM }>(i * PhaseSpace::PARTICLE_DIM) - } - - pub fn set_particle(&mut self, i: usize, position: Point, velocity: Vector) { - let mut view = self - .0 - .fixed_rows_mut::<{ PhaseSpace::PARTICLE_DIM }>(i * PhaseSpace::PARTICLE_DIM); - for i in 0..N { - view[i] = position[i]; - view[i + N] = velocity[i]; - } - } -} - -impl ParticleSystem { - fn collect_phase_space(&self) -> PhaseSpace { - let mut phase_space = PhaseSpace::new(self.particles.len()); - for (particle_index, particle) in self.particles.iter().enumerate() { - phase_space.set_particle(particle_index, particle.position, particle.velocity); - } - phase_space - } - - fn compute_derivative(&self) -> PhaseSpace { - let mut phase_space = PhaseSpace::new(self.particles.len()); - for (particle_index, particle) in self.particles.iter().enumerate() { - phase_space.set_particle( - particle_index, - particle.velocity.into(), - particle.force / particle.mass, - ); - } - phase_space - } - - fn scatter_phase_space(&mut self, phase_space: &PhaseSpace) { - for (particle_index, particle) in &mut self.particles.iter_mut().enumerate() { - let view = phase_space.particle_view(particle_index); - for i in 0..N { - particle.position[i] = view[i]; - particle.velocity[i] = view[i + N]; - } - } - } -} - -pub trait Solver { - fn step(&mut self, dt: Scalar); -} - -#[cfg(test)] -mod tests { - use super::{ParticleSystem, PhaseSpace, Point, Scalar, Solver, Vector}; - use crate::particle_system::Particle; - - #[test] - fn test_collect_phase_space() { - let system = ParticleSystem { - particles: vec![Particle::new(Point::new(2.0, 3.0), 1.0)], - constraints: vec![], - forces: vec![], - t: 0.0, - }; - let phase_space = system.collect_phase_space(); - - assert!( - !phase_space.0.is_empty(), - "Phase space has to contain non-zero values" - ); - } - - #[test] - fn test_scatter_phase_space() { - let mut phase_space = PhaseSpace::new(2); - phase_space.set_particle(1, Point::new(5.0, 7.0), Vector::x()); - - let mut system = ParticleSystem { - particles: vec![ - Particle::new(Point::new(0.0, 0.0), 1.0), - Particle::new(Point::new(0.0, 0.0), 1.0), - ], - constraints: vec![], - forces: vec![], - t: 0.0, - }; - - system.scatter_phase_space(&phase_space); - - assert!( - !system.particles[1].velocity.is_empty(), - "Velocity has to be set" - ); - assert!( - !system.particles[1].position.is_empty(), - "Position has to be set" - ); - } - - fn simulate_falling_ball( - fall_time: Scalar, - dt: Scalar, - mass: Scalar, - ) -> (Vector, Vector, ParticleSystem) { - let gravity = -9.8 * Vector::y(); - - let mut system = ParticleSystem { - particles: vec![Particle::new(Point::origin(), mass)], - constraints: vec![], - forces: vec![], - t: 0.0, - }; - - let iterations = (fall_time / dt) as usize; - - for _ in 0..iterations { - for particle in &mut system.particles { - particle.reset_force(); - particle.apply_force(gravity * particle.mass); - } - system.step(dt); - } - - let expected_velocity = gravity * fall_time; // vt - let expected_position = gravity * fall_time * fall_time / 2.0; // at^2 / 2 - - ( - system.particles[0].position.coords - expected_position, - system.particles[0].velocity - expected_velocity, - system, - ) - } - - #[test] - fn ball_should_fall() { - let (position_error, velocity_error, _) = simulate_falling_ball(10.0, 0.01, 3.0); - assert!( - position_error.norm() < 0.01, - "Position error is too high: {}", - position_error, - ); - assert!( - velocity_error.norm() < 0.01, - "Velocity error is too high: {}", - velocity_error, - ); - } - - #[test] - fn freefall_different_masses() { - let (_, _, system1) = simulate_falling_ball(10.0, 0.01, 2.0); - let (_, _, system2) = simulate_falling_ball(10.0, 0.01, 10.0); - - let diff = system1.particles[0].position - system2.particles[0].position; - assert!( - diff.norm() < 0.01, - "Points with different masses should fall with the same speed, diff: {}", - diff - ); - } -} -- cgit v1.2.3