summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Cargo.lock169
-rw-r--r--Cargo.toml2
-rw-r--r--README.md4
-rw-r--r--src/main.rs169
-rw-r--r--src/particle_system.rs204
5 files changed, 387 insertions, 161 deletions
diff --git a/Cargo.lock b/Cargo.lock
index 0971149..a26bd90 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -4,9 +4,9 @@ version = 3
[[package]]
name = "approx"
-version = "0.4.0"
+version = "0.5.1"
source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "3f2a05fd1bd10b2527e20a2cd32d8873d115b8b39fe219ee25f42a8aca6ba278"
+checksum = "cab112f0a86d568ea0e627cc1d6be74a1e9cd55214684db5561995f6dad897c6"
dependencies = [
"num-traits",
]
@@ -18,12 +18,84 @@ source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "ace50bade8e6234aa140d9a2f552bbee1db4d353f69b8217bc503490fc1a9f26"
[[package]]
-name = "cgmath"
-version = "0.18.0"
+name = "bytemuck"
+version = "1.20.0"
source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "1a98d30140e3296250832bbaaff83b27dcd6fa3cc70fb6f1f3e5c9c0023b5317"
+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",
]
@@ -37,8 +109,93 @@ dependencies = [
]
[[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 = [
- "cgmath",
+ "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
index 6f18a73..11c9611 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -4,4 +4,4 @@ version = "0.1.0"
edition = "2021"
[dependencies]
-cgmath = "0.18.0"
+nalgebra = "0.33.2"
diff --git a/README.md b/README.md
new file mode 100644
index 0000000..59d5c16
--- /dev/null
+++ b/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/src/main.rs b/src/main.rs
index ab2e15f..fcf19c7 100644
--- a/src/main.rs
+++ b/src/main.rs
@@ -1,167 +1,28 @@
-use std::io::Write;
-use std::{fs::File, path::PathBuf};
-use cgmath::{Angle, ElementWise, InnerSpace, Matrix3, MetricSpace, Rad, Vector3, Zero};
-
-type Vector = Vector3<f32>;
-
-#[derive(Debug, Clone, Copy)]
-struct Particle {
- position: Vector,
- position_old: Vector,
- // velocity: Vector,
- acceleration: Vector,
- mass: f32,
- fixed: bool,
-}
-
-impl Particle {
- pub fn new(position: Vector, mass: f32, fixed: bool) -> Self {
- Self {
- position,
- position_old: position,
- acceleration: Vector::zero(),
- mass,
- fixed,
- }
-
- }
-
- fn apply_force(&mut self, force: Vector) {
- self.acceleration += force / self.mass;
- }
-
- fn tick(&mut self, dt: f32) {
- if !self.fixed {
- // Verlet method
- let velocity = self.position - self.position_old;
- self.position_old = self.position;
- self.position += velocity + self.acceleration * dt * dt;
- }
-
- self.acceleration = Vector::zero();
- }
-}
-
-#[derive(Debug, Clone, Copy)]
-struct Spring {
- rest_length: f32,
- stiffness: f32,
- particle_ids: [usize; 2],
-}
-
-#[derive(Debug, Clone)]
-struct Simulation {
- particles: Vec<Particle>,
- springs: Vec<Spring>,
-}
-
-impl Simulation {
- fn tick(&mut self, dt: f32) {
- for spring in self.springs.iter() {
- let [a, b] = spring.particle_ids;
- let dist_vec = self.particles[b].position - self.particles[a].position;
- let force_magnitude = spring.stiffness * (dist_vec.magnitude() - spring.rest_length);
- let normalized = dist_vec.normalize();
-
- self.particles[a].apply_force(normalized.mul_element_wise(force_magnitude));
- self.particles[b].apply_force(normalized.mul_element_wise(force_magnitude * -1.0));
- }
- for p in &mut self.particles {
- p.tick(dt);
- }
- }
-}
-
-
-struct PPM {
- prefix: PathBuf,
- width: usize,
- height: usize,
- // pixels_per_unit: usize,
-}
-
-impl PPM {
- fn render_particles(&self, particles: &Vec<Particle>) -> 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 f32, 0.0, (pixel_row as f32) * -1.0) + Vector::new(self.width as f32 / -2.0, 0.0, self.height as f32 / 2.0);
- let color = match particles.iter().any(|p| p.position.distance(point) <= p.mass / 2.0) {
- true => black,
- false => white,
- };
- s += color;
- }
- }
- s
- }
-
- fn save_frame(&self, particles: &Vec<Particle>, time: f32) {
- 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();
- }
-}
+use particle_system::{Particle, ParticleSystem, Point, Vector};
+mod particle_system;
fn main() {
- let rotation_matrix = Matrix3::from_axis_angle(Vector::unit_y(), -Rad::full_turn() / 12.0);
-
- let mut simulation = Simulation {
- particles: vec![
- Particle::new(Vector::zero(), 5.0, true),
- Particle::new(rotation_matrix * Vector::unit_x() * 40.0, 30.0, false),
- Particle::new(rotation_matrix * Vector::unit_x() * 40.0 + rotation_matrix * rotation_matrix * Vector::unit_x() * 50.0, 15.0, false)
- ],
- springs: vec![
- Spring {
- particle_ids: [0, 1],
- rest_length: 40.0,
- stiffness: 400.0,
- },
- Spring {
- particle_ids: [1, 2],
- rest_length: 50.0,
- stiffness: 5.0,
- },
- ],
+ let dt = 0.01;
+ let mut system = ParticleSystem {
+ particles: vec![Particle::new(Point::origin(), 1.0)],
+ t: 0.0,
};
+ let gravity = -9.8 * Vector::y();
- let ppm = PPM {
- prefix: PathBuf::from("./out"),
- width: 150,
- height: 800,
- };
-
- let gravity = Vector::new(0.0, 0.0, -9.8);
-
- let time_step = 0.01;
- let mut tick = 1;
- let mut time = 0.0;
- while time < 130.0 {
- {
- let p = &mut simulation.particles[0];
- p.apply_force(gravity * -p.mass);
- }
- let mut all_particles = Vec::new();
- for p in &mut simulation.particles {
- p.apply_force(gravity * p.mass);
+ for i in 0..3 {
+ println!("Iteration #{i}");
+ for particle in &mut system.particles {
+ particle.reset_force();
}
- simulation.tick(time_step);
- all_particles.append(&mut simulation.particles.clone());
- if tick % 10 == 0 {
- ppm.save_frame(&all_particles, time);
+ for particle in &mut system.particles {
+ particle.apply_force(gravity);
}
+ system.euler_step(dt);
- time += time_step;
- tick +=1;
+ println!("{:?}", system);
}
}
diff --git a/src/particle_system.rs b/src/particle_system.rs
new file mode 100644
index 0000000..0556764
--- /dev/null
+++ b/src/particle_system.rs
@@ -0,0 +1,204 @@
+use nalgebra::{Const, DVector, Dyn, Matrix, Point as PointBase, SVector, ViewStorage};
+
+const N: usize = 2;
+pub type Scalar = f32;
+
+pub type Vector = SVector<Scalar, N>;
+pub type Point = PointBase<Scalar, N>;
+
+#[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()
+ }
+}
+
+/// A vector of concatenated position and velocity components of each particle
+#[derive(Debug)]
+pub struct PhaseSpace(DVector<Scalar>);
+type ParticleView<'a> = Matrix<
+ f32,
+ Const<{ PhaseSpace::PARTICLE_DIM }>,
+ Const<1>,
+ ViewStorage<'a, f32, 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::<Scalar>::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];
+ }
+ }
+}
+
+#[derive(Debug)]
+pub struct ParticleSystem {
+ pub particles: Vec<Particle>,
+
+ /// Simulation clock
+ pub t: Scalar,
+}
+
+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 fn euler_step(&mut self, dt: Scalar) {
+ let derivative = self.compute_derivative();
+ let mut state = self.collect_phase_space();
+
+ state.0 += derivative.0 * dt;
+ self.scatter_phase_space(&state);
+
+ self.t += dt;
+ }
+}
+
+#[cfg(test)]
+mod tests {
+ use super::{Particle, ParticleSystem, PhaseSpace, Point, Scalar, Vector};
+
+ #[test]
+ fn test_collect_phase_space() {
+ let system = ParticleSystem {
+ particles: vec![Particle::new(Point::new(2.0, 3.0), 1.0)],
+ 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),
+ ],
+ 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) -> (Vector, Vector) {
+ let gravity = -9.8 * Vector::y();
+
+ let mut system = ParticleSystem {
+ particles: vec![Particle::new(Point::origin(), 1.0)],
+ t: 0.0,
+ };
+
+ let iterations = (fall_time / dt) as usize;
+
+ for _ in 0..iterations {
+ for particle in &mut system.particles {
+ particle.reset_force();
+ }
+
+ for particle in &mut system.particles {
+ particle.apply_force(gravity);
+ }
+
+ system.euler_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,
+ )
+ }
+
+ #[test]
+ fn ball_should_fall() {
+ let (position_error, velocity_error) = simulate_falling_ball(10.0, 0.01);
+ assert!(position_error.norm() < 0.5, "Position error has is too high");
+ assert!(velocity_error.norm() < 0.5, "Velocity error has is too high");
+ }
+}
generated by cgit v1.2.3 (git 2.25.1) at 2025-05-10 14:17:07 +0000