Browse Source

Basic prototype

master
Stephen 5 months ago
parent
commit
d6d73c17aa
8 changed files with 346 additions and 91 deletions
  1. +15
    -0
      Cargo.lock
  2. +2
    -0
      Cargo.toml
  3. +20
    -0
      README.md
  4. +28
    -38
      src/barnes_hut.rs
  5. +97
    -44
      src/main.rs
  6. +6
    -9
      src/rigid_point.rs
  7. +104
    -0
      src/vector3d.rs
  8. +74
    -0
      src/video_3d.rs

+ 15
- 0
Cargo.lock View File

@ -83,6 +83,16 @@ source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "1b20b618342cf9891c292c4f5ac2cde7287cc5c87e87e9c769d617793607dec1"
[[package]]
name = "bincode"
version = "1.3.1"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "f30d3a39baa26f9651f17b375061f3233dde33424a8b72b0dbe93a68a0bc896d"
dependencies = [
"byteorder",
"serde",
]
[[package]]
name = "bitflags"
version = "1.2.1"
source = "registry+https://github.com/rust-lang/crates.io-index"
@ -311,10 +321,12 @@ checksum = "a06f77d526c1a601b7c4cdd98f54b5eaabffc14d5f2f0296febdc7f357c6d3ba"
name = "galaxy3d"
version = "0.1.0"
dependencies = [
"bincode",
"kiss3d",
"nalgebra",
"rand 0.7.3",
"rayon",
"serde",
]
[[package]]
@ -1228,6 +1240,9 @@ name = "serde"
version = "1.0.114"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "5317f7588f0a5078ee60ef675ef96735a1442132dc645eb1d12c018620ed8cd3"
dependencies = [
"serde_derive",
]
[[package]]
name = "serde_derive"


+ 2
- 0
Cargo.toml View File

@ -11,6 +11,8 @@ kiss3d = "0.24.1"
nalgebra = "0.21"
rand = "0.7"
rayon = "1.3"
serde = { version = "1.0", features = ["derive"] }
bincode = "1.3"
[target.x86_64-unknown-linux-gnu]
rustflags = [


+ 20
- 0
README.md View File

@ -0,0 +1,20 @@
# galaxy3d-rust
A simple 3D gravitational simulator I wrote in Rust.
## Usage
First, generate a 3D video.
`cargo run --release -- create`
This will generate a `video.3dv` file.
Next, watch the video.
`cargo run --release -- watch`
This will load the `video.3dv` file, so that you can watch it.
This is very much a WIP, and all parameters are hardcoded for now.

+ 28
- 38
src/barnes_hut.rs View File

@ -1,30 +1,29 @@
extern crate nalgebra as na;
use na::{Vector3};
use crate::rigid_point::RigidPoint;
use crate::vector3d::Vector3D;
const G: f64 = 0.00667408;
const MAX_DEPTH: i32 = 100;
const THETA: f64 = 0.5;
pub struct Octree {
nodes: Vec<OctNode>
}
pub struct OctNode {
parent: isize, // -1 = no parent
children: [isize; 8],
parent: isize, // -1 = no parent. TODO get rid of all the parent stuff - we never use it
children: [isize; 8], // TODO make this an Option<[usize; 8]>
mass: f64,
min_coord: Vector3<f64>,
max_coord: Vector3<f64>,
min_coord: Vector3D,
max_coord: Vector3D,
num_particles: usize,
avg_weighted_position: Vector3<f64>,// Divide by mass for this to be accurate!
center_of_mass: Vector3<f64> // Calculated once from avg_weighted_position
avg_weighted_position: Vector3D, // Divide by mass for this to be accurate!
center_of_mass: Vector3D // Calculated once from avg_weighted_position
}
impl Octree {
pub fn new(min: Vector3<f64>, max: Vector3<f64>) -> Self {
pub fn new(min: Vector3D, max: Vector3D) -> Self {
let mut nodes = Vec::new();
nodes.push(OctNode::new(min, max, -1)); // Create our root node
Octree { nodes: nodes }
@ -51,24 +50,27 @@ impl Octree {
let min_coord = self.nodes[cur_node].min_coord;
let max_coord = self.nodes[cur_node].max_coord;
let width = max_coord.x - min_coord.x;
let height = max_coord.y - min_coord.y;
let depth = max_coord.z - min_coord.z;
let s = if width > height { width } else { height };
let s = if depth > s { depth } else { s };
let s_squared = s * s;
//let s_squared = (self.nodes[cur_node].max_coord -
// self.nodes[cur_node].min_coord).norm_squared(); // width^2
let dis_vec = self.nodes[cur_node].center_of_mass - particle.position;
let d_squared = dis_vec.norm_squared();
let theta_squared = 0.5 * 0.5; // TODO make this const
let d_squared = dis_vec.dist_squared();
let theta_squared = THETA * THETA;
let ratio = s_squared / d_squared;
if self.nodes[cur_node].num_particles == 1 || ratio < theta_squared {
println!("s/d = {}", ratio.sqrt());
//We can approximate
if(d_squared < 0.1) { return; }
if d_squared < 0.1 { return; } // Worse, but higher-performance, method of smoothing
let other = &self.nodes[cur_node];
let dist = d_squared.sqrt();
let v = G * other.mass * delta / d_squared; // / (d_squared + 0.1).sqrt();
let v = G * other.mass * delta / d_squared; // / (d_squared + 0.1).sqrt(); // Uncomment this to enable smoothing for short distances
particle.velocity += v * dis_vec / dist;
}
else {
@ -113,9 +115,6 @@ impl Octree {
let mut particle_octs: [Vec<RigidPoint>; 8] = [vec![], vec![], vec![], vec![],
vec![], vec![], vec![], vec![]];
for particle in particles {
/*println!("// {} {} {} //", particle.position, center_point, self.nodes[cur_node].max_coord);
let mut input = String::new();
io::stdin().read_line(&mut input);*/
// Update ourself
self.nodes[cur_node].mass += particle.mass;
@ -134,7 +133,7 @@ impl Octree {
}
// TODO wire this up
pub fn get_id_from_center(center: Vector3<f64>, point: Vector3<f64>) -> usize {
pub fn get_id_from_center(center: Vector3D, point: Vector3D) -> usize {
let offset = point - center;
//We can look at the sign of the components of the offset to figure out which octant it is in.
let x_offset = if offset.x > 0.0 { 0 } else { 1 };
@ -143,12 +142,12 @@ impl Octree {
x_offset * 1 + y_offset * 2 + z_offset * 4 // basic binary stuff here
}
pub fn get_octant_bounding_box_from_id(min: Vector3<f64>,
max: Vector3<f64>,
center: Vector3<f64>,
idx: usize) -> (Vector3<f64>, Vector3<f64>) {
let mut min_coord: Vector3<f64> = Vector3::new(0.0, 0.0, 0.0);
let mut max_coord: Vector3<f64> = Vector3::new(0.0, 0.0, 0.0);
pub fn get_octant_bounding_box_from_id(min: Vector3D,
max: Vector3D,
center: Vector3D,
idx: usize) -> (Vector3D, Vector3D) {
let mut min_coord: Vector3D = Vector3D::new(0.0, 0.0, 0.0);
let mut max_coord: Vector3D = Vector3D::new(0.0, 0.0, 0.0);
if (idx & 1) != 0 {
min_coord.x = min.x;
max_coord.x = center.x;
@ -178,19 +177,10 @@ impl Octree {
(min_coord, max_coord)
}
pub fn check_for_collision(&self, particle: RigidPoint) {
check_for_collision_with_specific_node(particle, 0);
}
fn check_for_collision_with_specific_node(&self, particle: RigidPoint, idx: usize) {
let cur_node = self.nodes[idx];
}
}
impl OctNode {
fn new(min: Vector3<f64>, max: Vector3<f64>, parent: isize) -> Self {
fn new(min: Vector3D, max: Vector3D, parent: isize) -> Self {
OctNode {
parent: parent,
children: [-1, -1, -1, -1, -1, -1, -1, -1],
@ -198,8 +188,8 @@ impl OctNode {
min_coord: min,
max_coord: max,
num_particles: 0,
avg_weighted_position: Vector3::new(0.0, 0.0, 0.0),
center_of_mass: Vector3::new(0.0, 0.0, 0.0)
avg_weighted_position: Vector3D::new(0.0, 0.0, 0.0),
center_of_mass: Vector3D::new(0.0, 0.0, 0.0)
}
}
}

+ 97
- 44
src/main.rs View File

@ -1,21 +1,23 @@
extern crate kiss3d;
extern crate nalgebra as na;
use na::{Vector3, Translation3, Point3, Point2};
use na::{Vector3, Translation3, Point3};
use kiss3d::window::Window;
use kiss3d::light::Light;
use kiss3d::scene::SceneNode;
use kiss3d::text::Font;
use rand::prelude::*;
use std::time::Instant;
use rayon::prelude::*;
use std::sync::{Arc, Mutex};
use std::env;
mod barnes_hut;
mod rigid_point;
mod video_3d;
mod vector3d;
use rigid_point::RigidPoint;
use barnes_hut::Octree;
use vector3d::Vector3D;
use video_3d::*;
const G: f64 = 0.00667408;
@ -27,19 +29,17 @@ struct Line {
struct Universe {
particles: Vec<RigidPoint>,
nodes: Vec<SceneNode>,
lines: Vec<Line>,
delta: f64,
}
impl Universe {
fn new(delta: f64) -> Self {
Universe { particles: Vec::new(), nodes: Vec::new(), lines: Vec::new(), delta: delta }
Universe { particles: Vec::new(), lines: Vec::new(), delta: delta }
}
fn push(&mut self, particle: RigidPoint, node: SceneNode) {
fn push(&mut self, particle: RigidPoint) {
self.particles.push(particle);
self.nodes.push(node)
}
fn add_line(&mut self, start: Point3<f32>, end: Point3<f32>) {
@ -57,11 +57,11 @@ impl Universe {
let mut energy: f64 = 0.0;
for a in self.particles.clone() {
//Calculate kinetic energy
energy += 0.5 * a.mass * a.velocity.norm_squared();
energy += 0.5 * a.mass * a.velocity.dist_squared();
//Calculate gravitational potential energy
for b in &self.particles {
if a.index != b.index {
let dist = na::distance(&Point3::from(a.position), &Point3::from(b.position));
let dist = (a.position - b.position).dist();
let energy_one = -G * b.mass * a.mass / dist;
// println!("dist: {} energy: {}", dist, energy_one);
energy += energy_one;
@ -92,15 +92,15 @@ impl Universe {
let new_r = (r1 * r1 * r1 + r2 * r2 * r2).cbrt();
//Destroy i
self.nodes[i].unlink();
//self.nodes[i].unlink();
self.particles[i].index = -1;
//Update j
self.nodes[j].unlink();
/*self.nodes[j].unlink();
let mut new_node = window.add_sphere(new_r as f32);
new_node.set_color(1.0, 0.0, 0.0);
new_node.append_translation(&Translation3::from(b.position_f32()));
self.nodes[j] = new_node;
self.nodes[j] = new_node;*/
self.particles[j].position = new_pos;
self.particles[j].velocity = new_vel;
self.particles[j].radius = new_r;
@ -110,44 +110,92 @@ impl Universe {
}
fn main() {
let mut window = Window::new("Galaxy simulator");
let args: Vec<String> = env::args().collect();
if args.len() < 2 {
println!("Usage: ./{} watch|create", args[0]);
return;
}
match args[1].as_ref() {
"watch" => watch_video(),
"create" => create_video(),
_ => {
println!("Usage: ./{} watch|create", args[0]);
return;
}
}
}
window.set_light(Light::StickToCamera);
fn watch_video() {
let vid = Video3D::load_file("video.3dv");
let mut window = Window::new("3DVideo viewer");
let mut i: u64 = 0;
let mut nodes: Vec<SceneNode> = Vec::new();
while window.render() {
vid.get_frame(i as usize % vid.frames.len()).spheres.into_iter().enumerate().for_each(|(j, sphere)| {
if j >= nodes.len() {
let mut node = window.add_sphere(sphere.radius as f32);
node.set_color(sphere.r as f32, sphere.g as f32, sphere.b as f32);
node.append_translation(&Translation3::new(sphere.x as f32, sphere.y as f32, sphere.z as f32));
nodes.push(node);
}
else {
// We can reuse existing spheres
nodes[j].set_color(sphere.r as f32, sphere.g as f32, sphere.b as f32);
nodes[j].set_local_translation(Translation3::new(sphere.x as f32, sphere.y as f32, sphere.z as f32));
}
});
let delta = 0.005;
let mut universe = generate_random_points(&mut window, delta);
// TODO destroy nodes that aren't being reused
i += 1;
}
}
//let mut start_time = Instant::now();
while window.render() {
fn create_video() {
/* ====== RENDER TO A 3D VIDEO ====== */
let mut vid = Video3D::new(60);
let delta = 0.005;
let mut universe = generate_random_points(delta);
for iteration in 0..100 {
println!("Start iter");
for _ in 0..10 {
tick(&mut universe, &mut window);
tick(&mut universe);
calc_velocities(&mut universe, delta);
}
universe.render(&mut window);
//Calculate energy. Used to tell how much the algorithm sucks
/*window.draw_text(&format!("Total energy: {}", universe.calc_total_energy()),
&Point2::new(0.0, 0.0),
120.0,
&Font::default(),
&Point3::new(1.0, 1.0, 1.0));*/
//start_time = Instant::now();
}
println!("Stepping complete");
//Create frame from universe
let mut frame = Frame3D::new();
for p in &universe.particles {
let pos = p.position;
frame.push_sphere(Sphere3D::new(pos.x, pos.y, pos.z,
p.radius,
1.0, 1.0, 1.0));
}
//Push frame to vid
vid.push_frame(frame);
//if iteration % 10 == 0 {
println!("{}% complete", iteration / 10);
//}
}
vid.save_file("video.3dv");
}
fn generate_random_points(window: &mut Window, delta: f64) -> Universe {
fn generate_random_points(delta: f64) -> Universe {
let mut ret: Universe = Universe::new(delta);
let mut rng = rand::thread_rng();
for i in 0..5_000 {
for i in 0..100_000 {
let x: f64 = rng.gen_range(-50.0, 50.0);
let y: f64 = rng.gen_range(-50.0, 50.0);
let z: f64 = rng.gen_range(-5.0, 5.0); // Flat galaxy
let position: Vector3<f64> = Vector3::new(x, y, z);
let velocity: Vector3<f64> = Vector3::new(-y, x, 0.0);
let position: Vector3D = Vector3D::new(x, y, z);
let velocity: Vector3D = Vector3D::new(-y, x, 0.0);
// let velocity: Vector3<f64> = Vector3::new(0.0, 0.0, 0.0);
let node = window.add_sphere(0.1);
ret.push(RigidPoint::new(position, velocity / 60.0, 80.0, i, 0.1), node);
ret.push(RigidPoint::new(position, velocity / 60.0, 0.1, i, 0.1));
}
//Move velocities forward slightly
@ -160,7 +208,7 @@ fn generate_random_points(window: &mut Window, delta: f64) -> Universe {
fn calc_velocities(u: &mut Universe, delta: f64) {
//Find min/max
//println!("Finding min/max");
println!("Finding min/max");
let mut min: Vector3<f64> = Vector3::new(0.0, 0.0, 0.0);
let mut max: Vector3<f64> = Vector3::new(0.0, 0.0, 0.0);
for particle in &u.particles {
@ -189,18 +237,19 @@ fn calc_velocities(u: &mut Universe, delta: f64) {
}
}
// println!("Constructing tree");
println!("Constructing tree");
let mut oct = Octree::new(min, max);
let mut oct = Octree::new(Vector3D::from_na_vector(min),
Vector3D::from_na_vector(max));
oct.construct_tree(u.particles.clone().into_iter().filter(|x| {x.index != -1}).collect());
//println!("Barnes hut");
println!("Barnes hut");
u.particles.par_iter_mut().for_each(|mut p| {
oct.barnes_hut(&mut p, delta);
});
//println!("Done");
println!("Done");
}
/*fn calc_velocities(u: &mut Universe, delta: f64) {
@ -216,6 +265,7 @@ fn calc_velocities(u: &mut Universe, delta: f64) {
});
}*/
/*
fn calc_gravity(body: &mut RigidPoint, other: &RigidPoint, delta: f64) {
if body.index == other.index || body.index == -1 || other.index == -1 { // can't be the same
return;
@ -231,10 +281,11 @@ fn calc_gravity(body: &mut RigidPoint, other: &RigidPoint, delta: f64) {
// let v = -G * other.mass / dist_squared * delta;
let v = -G * other.mass * delta / dist / (dist_squared + 0.1).sqrt(); // softening
let hyp = dist;
body.velocity += Vector3::new(v * dx / hyp, v * dy / hyp, v * dz / hyp);
body.velocity += Vector3D::new(v * dx / hyp, v * dy / hyp, v * dz / hyp);
}
*/
fn tick(u: &mut Universe, window: &mut Window) {
fn tick(u: &mut Universe) {
let delta = u.delta;
u.particles.par_iter_mut().for_each(|a| {
if a.index == -1 {
@ -243,11 +294,13 @@ fn tick(u: &mut Universe, window: &mut Window) {
a.position += a.velocity * delta;
});
/*
u.particles.clone().iter().for_each(|a| {
if a.index != -1 {
a.update_node(&mut u.nodes[a.index as usize]);
}
});
*/
// Add line for particle 0
/*{


+ 6
- 9
src/rigid_point.rs View File

@ -1,19 +1,20 @@
extern crate nalgebra as na;
use na::{Vector3, Translation3, Point3, Point2};
use kiss3d::scene::SceneNode;
use na::Vector3;
use crate::vector3d::Vector3D;
#[derive(Clone)]
pub struct RigidPoint {
pub position: Vector3<f64>,
pub velocity: Vector3<f64>, // Velocity is one half-timestep in the future
pub position: Vector3D,
pub velocity: Vector3D, // Velocity is one half-timestep in the future
pub mass: f64,
pub index: isize,
pub radius: f64
}
impl RigidPoint {
pub fn new(position: Vector3<f64>, velocity: Vector3<f64>, mass: f64, index: isize, radius: f64) -> Self {
pub fn new(position: Vector3D,
velocity: Vector3D, mass: f64, index: isize, radius: f64) -> Self {
RigidPoint {
position: position,
velocity: velocity,
@ -23,10 +24,6 @@ impl RigidPoint {
}
}
pub fn update_node(&self, node: &mut SceneNode) {
node.set_local_translation(Translation3::from(self.position_f32()));
}
pub fn position_f32(&self) -> Vector3<f32> {
Vector3::new(self.position.x as f32, self.position.y as f32, self.position.z as f32)
}


+ 104
- 0
src/vector3d.rs View File

@ -0,0 +1,104 @@
use std::ops::{
Add,
Sub,
Mul,
Div,
AddAssign
};
use na::Vector3;
// Fast vectors
#[derive(Clone, Copy)]
pub struct Vector3D {
pub x: f64,
pub y: f64,
pub z: f64
}
impl Vector3D {
pub fn new(x: f64, y: f64, z: f64) -> Self {
Self { x: x, y: y, z: z }
}
fn get_na_vector(&self) -> Vector3<f64> {
Vector3::new(self.x, self.y, self.z)
}
pub fn from_na_vector(v: Vector3<f64>) -> Self {
Self {
x: v.x as f64,
y: v.y as f64,
z: v.z as f64
}
}
pub fn dist_squared(self) -> f64 {
self.x * self.x +
self.y * self.y +
self.z * self.z
}
pub fn dist(self) -> f64 {
self.dist_squared().sqrt()
}
}
impl Mul<f64> for Vector3D {
type Output = Self;
fn mul(self, rhs: f64) -> Self {
Self { x: self.x * rhs, y: self.y * rhs, z: self.z * rhs }
}
}
impl Mul<Vector3D> for f64 {
type Output = Vector3D;
fn mul(self, rhs: Vector3D) -> Vector3D {
Vector3D {
x: self * rhs.x,
y: self * rhs.y,
z: self * rhs.z
}
}
}
impl Div<f64> for Vector3D {
type Output = Self;
fn div(self, rhs: f64) -> Self {
Self { x: self.x / rhs, y: self.y / rhs, z: self.z / rhs }
}
}
impl Sub for Vector3D {
type Output = Self;
fn sub(self, rhs: Self) -> Self {
Self {
x: self.x - rhs.x,
y: self.y - rhs.y,
z: self.z - rhs.z
}
}
}
impl Add for Vector3D {
type Output = Self;
fn add(self, other: Self) -> Self {
Self {
x: self.x + other.x,
y: self.y + other.y,
z: self.z + other.z
}
}
}
impl AddAssign for Vector3D {
fn add_assign(&mut self, other: Self) {
self.x += other.x;
self.y += other.y;
self.z += other.z;
}
}

+ 74
- 0
src/video_3d.rs View File

@ -0,0 +1,74 @@
use serde::{Serialize, Deserialize};
use std::io::prelude::*;
use std::io::BufWriter;
use std::fs::File;
#[derive(Serialize, Deserialize, Debug)]
pub struct Video3D {
framerate: u32,
pub frames: Vec<Frame3D>
}
impl Video3D {
pub fn new(framerate: u32) -> Self {
Self { framerate, frames: Vec::new() }
}
pub fn push_frame(&mut self, f: Frame3D) {
self.frames.push(f);
}
pub fn save_file(&self, filename: &str) -> std::io::Result<()> {
let mut buffer = BufWriter::new(File::create(filename)?);
buffer.write_all(&bincode::serialize(self).unwrap())?;
buffer.flush()?;
Ok(())
}
pub fn load_file(filename: &str) -> Self {
let mut buffer = Vec::new();
File::open(filename).unwrap().read_to_end(&mut buffer).unwrap();
bincode::deserialize(&mut buffer).unwrap()
}
pub fn get_frame(&self, i: usize) -> Frame3D {
self.frames[i].clone()
}
}
// May add other types in the future
#[derive(Serialize, Deserialize, Debug, Clone)]
pub struct Frame3D {
pub spheres: Vec<Sphere3D>
}
impl Frame3D {
pub fn new() -> Self {
Self { spheres: Vec::new() }
}
pub fn push_sphere(&mut self, s: Sphere3D) {
self.spheres.push(s);
}
}
#[derive(Serialize, Deserialize, Debug, Clone)]
pub struct Sphere3D {
pub x: f64,
pub y: f64,
pub z: f64,
pub radius: f64,
// Range from 0.0 to 1.0
pub r: f64,
pub g: f64,
pub b: f64,
}
impl Sphere3D {
// TODO make this signature less ridiculous
pub fn new(x: f64, y: f64, z: f64, radius: f64, r: f64, g: f64, b: f64) -> Self {
Self { x, y, z, radius, r, g, b }
}
}

Loading…
Cancel
Save