Björn Steinbrink 8 months ago
parent
commit
5dc0af89ca
3 changed files with 49 additions and 26 deletions
  1. +4
    -1
      Cargo.toml
  2. +42
    -23
      src/barnes_hut.rs
  3. +3
    -2
      src/main.rs

+ 4
- 1
Cargo.toml View File

@ -17,4 +17,7 @@ bincode = "1.3"
[target.x86_64-unknown-linux-gnu]
rustflags = [
"-C", "link-arg=-fuse-ld=lld",
]
]
[profile.release]
debug = 1

+ 42
- 23
src/barnes_hut.rs View File

@ -18,27 +18,35 @@ pub struct Octree {
}
impl Octree {
pub fn new(min: Vector3D, max: Vector3D) -> Self {
pub fn new(min: Vector3D, max: Vector3D, n: usize) -> Self {
let mut o = Octree {
num_particles: Vec::with_capacity(n),
children: Vec::with_capacity(n),
mass: Vec::with_capacity(n),
min_coord: Vec::with_capacity(n),
max_coord: Vec::with_capacity(n),
avg_weighted_position: Vec::with_capacity(n),
center_of_mass: Vec::with_capacity(n),
s: Vec::with_capacity(n),
};
let width = max.x - min.x;
let height = max.y - min.y;
let depth = max.z - min.z;
let s = if width > height { width } else { height };
let s = if depth > s { depth } else { s };
Octree {
num_particles: vec![0],
children: vec![[-1; 8]],
mass: vec![0.0],
min_coord: vec![min],
max_coord: vec![max],
avg_weighted_position: vec![Vector3D::new(0.0, 0.0, 0.0)],
center_of_mass: Vec::new(),
s: vec![s],
}
o.num_particles.push(0);
o.children.push([-1; 8]);
o.mass.push(0.0);
o.min_coord.push(min);
o.max_coord.push(max);
o.avg_weighted_position.push(Vector3D::new(0.0, 0.0, 0.0));
o.s.push(s);
o
}
pub fn construct_tree(&mut self, particles: Vec<RigidPoint>) {
self.add_particles_to_specific_node(particles, 0);
pub fn construct_tree(&mut self, particles: &[&RigidPoint]) {
self.add_particles_to_specific_node(&particles, 0);
//calculate center of masses
for (&mass, &awp) in self.mass.iter().zip(&self.avg_weighted_position) {
@ -56,14 +64,17 @@ impl Octree {
delta: f64) {
if self.num_particles[cur_node] == 0 { return; }
let s = self.s[cur_node];
let s_squared = s * s;
let dis_vec = self.center_of_mass[cur_node] - particle.position;
let d_squared = dis_vec.dist_squared();
let theta_squared = THETA * THETA;
let ratio = s_squared / d_squared;
if self.num_particles[cur_node] == 1 || ratio < theta_squared {
if self.num_particles[cur_node] == 1 || {
let s = self.s[cur_node];
let s_squared = s * s;
let ratio = s_squared / d_squared;
let theta_squared = THETA * THETA;
ratio < theta_squared
}
{
// println!("s/d = {}", ratio.sqrt());
//We can approximate
if d_squared < 0.1 { return; } // Worse, but higher-performance, method of smoothing
@ -81,7 +92,7 @@ impl Octree {
}
}
fn add_particles_to_specific_node(&mut self, particles: Vec<RigidPoint>, cur_node: usize) {
fn add_particles_to_specific_node(&mut self, particles: &[&RigidPoint], cur_node: usize) {
if particles.len() == 0 {
return;
}
@ -123,9 +134,17 @@ impl Octree {
println!("This should never happen :S");
}
let mut particle_octs: [Vec<RigidPoint>; 8] = [vec![], vec![], vec![], vec![],
vec![], vec![], vec![], vec![]];
for particle in particles {
let mut particle_octs: [Vec<&RigidPoint>; 8] = [
Vec::with_capacity(particles.len()/4),
Vec::with_capacity(particles.len()/4),
Vec::with_capacity(particles.len()/4),
Vec::with_capacity(particles.len()/4),
Vec::with_capacity(particles.len()/4),
Vec::with_capacity(particles.len()/4),
Vec::with_capacity(particles.len()/4),
Vec::with_capacity(particles.len()/4),
];
for &particle in particles {
// Update ourself
self.mass[cur_node] += particle.mass;
@ -138,7 +157,7 @@ impl Octree {
//Recurse
for (i, particle_oct) in particle_octs.iter().enumerate() {
self.add_particles_to_specific_node(particle_oct.clone(),
self.add_particles_to_specific_node(&particle_oct,
self.children[cur_node][i] as usize);
}
}


+ 3
- 2
src/main.rs View File

@ -239,9 +239,10 @@ fn calc_velocities(u: &mut Universe, delta: f64) {
println!("Constructing tree");
let p = u.particles.iter().filter(|x| {x.index != -1}).collect::<Vec<_>>();
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());
Vector3D::from_na_vector(max), p.len());
oct.construct_tree(&p);
println!("Barnes hut");


Loading…
Cancel
Save