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 --- 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 ----------------------------------------------- 14 files changed, 753 deletions(-) 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 (limited to 'src') 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