summaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/algebra/mod.rs1
-rw-r--r--src/algebra/subspace.rs69
-rw-r--r--src/constraint/anchor.rs42
-rw-r--r--src/constraint/beam.rs46
-rw-r--r--src/constraint/mod.rs139
-rw-r--r--src/constraint/slider.rs50
-rw-r--r--src/force/drag.rs15
-rw-r--r--src/force/gravity.rs15
-rw-r--r--src/force/mod.rs24
-rw-r--r--src/main.rs59
-rw-r--r--src/particle_system.rs47
-rw-r--r--src/ppm.rs41
-rw-r--r--src/solver/midpoint.rs18
-rw-r--r--src/solver/mod.rs187
14 files changed, 0 insertions, 753 deletions
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<Scalar, N, N>;
-
-pub struct Subspace<const DIM: usize> {
- point: Point,
- vectors: [Vector; DIM],
- projection_matrix: ProjectionMatrix,
-}
-
-pub type Line = Subspace<1>;
-pub type Plane = Subspace<2>;
-
-impl<const DIM: usize> Subspace<DIM> {
- 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::<Scalar, N, DIM>::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<Scalar>,
-}
-
-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<usize> {
- vec![self.particle_id]
- }
-
- fn c(&self, q: &DVector<Scalar>) -> Scalar {
- let position = q.fixed_rows(self.particle_id * N);
- (position - self.anchor.coords).norm()
- }
-
- fn set_jacobian(&mut self, jacobian: RowDVector<Scalar>) {
- self.jacobian = jacobian
- }
-
- fn jacobian_prev(&self) -> RowDVector<Scalar> {
- 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<Scalar>,
-}
-
-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<usize> {
- Vec::from(self.particle_ids)
- }
-
- fn c(&self, q: &DVector<Scalar>) -> Scalar {
- let a = q.fixed_rows::<N>(self.particle_ids[0] * N);
- let b = q.fixed_rows::<N>(self.particle_ids[1] * N);
-
- (a - b).norm() - self.length
- }
-
- fn set_jacobian(&mut self, jacobian: RowDVector<Scalar>) {
- self.jacobian = jacobian
- }
-
- fn jacobian_prev(&self) -> RowDVector<Scalar> {
- 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<usize>;
-
- fn c(&self, q: &DVector<Scalar>) -> Scalar;
-
- fn jacobian_prev(&self) -> RowDVector<Scalar>;
- fn set_jacobian(&mut self, jacobian: RowDVector<Scalar>);
-
- fn partial_derivative(&self, q: &DVector<Scalar>) -> RowDVector<Scalar> {
- 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<Scalar>,
- q_prev: &DVector<Scalar>,
- dt: Scalar,
- ) -> (Scalar, Scalar, RowDVector<Scalar>, RowDVector<Scalar>) {
- 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::<Scalar>::zeros(self.constraints.len(), size);
- let mut jacobian_dot = DMatrix::<Scalar>::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<Scalar>,
-}
-
-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<usize> {
- vec![self.particle_id]
- }
-
- fn c(&self, q: &DVector<Scalar>) -> Scalar {
- let position = q.fixed_rows::<N>(self.particle_id * N);
- let point = Point::from_coordinates(position.into());
- self.line.distance(point)
- }
-
- fn set_jacobian(&mut self, jacobian: RowDVector<Scalar>) {
- self.jacobian = jacobian
- }
-
- fn jacobian_prev(&self) -> RowDVector<Scalar> {
- 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<crate::particle_system::Particle>) {
- 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<crate::particle_system::Particle>) {
- 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<Particle>);
-}
-
-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<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()
- }
-}
-
-// #[derive(Debug)]
-pub struct ParticleSystem {
- pub particles: Vec<Particle>,
- pub constraints: Vec<Box<dyn Constraint>>,
- pub forces: Vec<Box<dyn Force>>,
-
- /// 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<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 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<Particle>, 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<Scalar>);
-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::<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];
- }
- }
-}
-
-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
- );
- }
-}
generated by cgit v1.2.3 (git 2.25.1) at 2025-05-10 12:34:09 +0000