|
|
|
@ -7,34 +7,42 @@ const G: f64 = 0.00667408; |
|
|
|
const THETA: f64 = 0.5;
|
|
|
|
|
|
|
|
pub struct Octree {
|
|
|
|
nodes: Vec<OctNode>
|
|
|
|
}
|
|
|
|
|
|
|
|
pub struct OctNode {
|
|
|
|
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: Vector3D,
|
|
|
|
max_coord: Vector3D,
|
|
|
|
num_particles: usize,
|
|
|
|
avg_weighted_position: Vector3D, // Divide by mass for this to be accurate!
|
|
|
|
center_of_mass: Vector3D // Calculated once from avg_weighted_position
|
|
|
|
num_particles: Vec<usize>,
|
|
|
|
children: Vec<[isize; 8]>, // TODO make this an Option<[usize; 8]>
|
|
|
|
mass: Vec<f64>,
|
|
|
|
min_coord: Vec<Vector3D>,
|
|
|
|
max_coord: Vec<Vector3D>,
|
|
|
|
avg_weighted_position: Vec<Vector3D>, // Divide by mass for this to be accurate!
|
|
|
|
center_of_mass: Vec<Vector3D>, // Calculated once from avg_weighted_position
|
|
|
|
s: Vec<f64>,
|
|
|
|
}
|
|
|
|
|
|
|
|
impl Octree {
|
|
|
|
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 }
|
|
|
|
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],
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
pub fn construct_tree(&mut self, particles: Vec<RigidPoint>) {
|
|
|
|
self.add_particles_to_specific_node(particles, 0);
|
|
|
|
|
|
|
|
//calculate center of masses
|
|
|
|
for n in &mut self.nodes {
|
|
|
|
n.center_of_mass = n.avg_weighted_position / n.mass;
|
|
|
|
for (&mass, &awp) in self.mass.iter().zip(&self.avg_weighted_position) {
|
|
|
|
self.center_of_mass.push(awp / mass);
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
@ -46,39 +54,29 @@ impl Octree { |
|
|
|
particle: &mut RigidPoint,
|
|
|
|
cur_node: usize,
|
|
|
|
delta: f64) {
|
|
|
|
let c_node = &self.nodes[cur_node];
|
|
|
|
if c_node.num_particles == 0 { return; }
|
|
|
|
|
|
|
|
let min_coord = c_node.min_coord;
|
|
|
|
let max_coord = c_node.max_coord;
|
|
|
|
if self.num_particles[cur_node] == 0 { return; }
|
|
|
|
|
|
|
|
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 = self.s[cur_node];
|
|
|
|
let s_squared = s * s;
|
|
|
|
let dis_vec = self.nodes[cur_node].center_of_mass - particle.position;
|
|
|
|
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 c_node.num_particles == 1 || ratio < theta_squared {
|
|
|
|
if self.num_particles[cur_node] == 1 || ratio < theta_squared {
|
|
|
|
// println!("s/d = {}", ratio.sqrt());
|
|
|
|
//We can approximate
|
|
|
|
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(); // Uncomment this to enable smoothing for short distances
|
|
|
|
let v = G * self.mass[cur_node] * delta / d_squared; // / (d_squared + 0.1).sqrt(); // Uncomment this to enable smoothing for short distances
|
|
|
|
particle.velocity += v * dis_vec / dist;
|
|
|
|
}
|
|
|
|
else {
|
|
|
|
//println!("A{}", self.nodes[cur_node].num_particles);
|
|
|
|
//Divide and conquer
|
|
|
|
for i in 0..8 {
|
|
|
|
self.barnes_hut_specific_node(particle, c_node.children[i] as usize, delta);
|
|
|
|
for &child in &self.children[cur_node] {
|
|
|
|
self.barnes_hut_specific_node(particle, child as usize, delta);
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
@ -89,24 +87,36 @@ impl Octree { |
|
|
|
}
|
|
|
|
|
|
|
|
if particles.len() == 1 {
|
|
|
|
self.nodes[cur_node].mass += particles[0].mass;
|
|
|
|
self.nodes[cur_node].num_particles += 1;
|
|
|
|
self.nodes[cur_node].avg_weighted_position += particles[0].position * particles[0].mass;
|
|
|
|
self.mass[cur_node] += particles[0].mass;
|
|
|
|
self.num_particles[cur_node] += 1;
|
|
|
|
self.avg_weighted_position[cur_node] += particles[0].position * particles[0].mass;
|
|
|
|
return;
|
|
|
|
}
|
|
|
|
|
|
|
|
let center_point = (self.nodes[cur_node].min_coord + self.nodes[cur_node].max_coord) / 2.0;
|
|
|
|
let center_point = (self.min_coord[cur_node] + self.max_coord[cur_node]) / 2.0;
|
|
|
|
|
|
|
|
if self.nodes[cur_node].children[0] == -1 {
|
|
|
|
if self.children[cur_node][0] == -1 {
|
|
|
|
for i in 0..8 {
|
|
|
|
let idx = self.nodes.len();
|
|
|
|
let (min, max) = Octree::get_octant_bounding_box_from_id(self.nodes[cur_node].min_coord,
|
|
|
|
self.nodes[cur_node].max_coord,
|
|
|
|
let idx = self.num_particles.len();
|
|
|
|
let (min, max) = Octree::get_octant_bounding_box_from_id(self.min_coord[cur_node],
|
|
|
|
self.max_coord[cur_node],
|
|
|
|
center_point,
|
|
|
|
i);
|
|
|
|
|
|
|
|
self.nodes.push(OctNode::new(min, max, cur_node as isize));
|
|
|
|
self.nodes[cur_node].children[i] = idx as isize;
|
|
|
|
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 };
|
|
|
|
|
|
|
|
self.num_particles.push(0);
|
|
|
|
self.children.push([-1; 8]);
|
|
|
|
self.mass.push(0.0);
|
|
|
|
self.min_coord.push(min);
|
|
|
|
self.max_coord.push(max);
|
|
|
|
self.avg_weighted_position.push(Vector3D::new(0.0, 0.0, 0.0));
|
|
|
|
self.s.push(s);
|
|
|
|
self.children[cur_node][i] = idx as isize;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
else {
|
|
|
|
@ -118,9 +128,9 @@ impl Octree { |
|
|
|
for particle in particles {
|
|
|
|
|
|
|
|
// Update ourself
|
|
|
|
self.nodes[cur_node].mass += particle.mass;
|
|
|
|
self.nodes[cur_node].num_particles += 1;
|
|
|
|
self.nodes[cur_node].avg_weighted_position += particle.position * particle.mass;
|
|
|
|
self.mass[cur_node] += particle.mass;
|
|
|
|
self.num_particles[cur_node] += 1;
|
|
|
|
self.avg_weighted_position[cur_node] += particle.position * particle.mass;
|
|
|
|
|
|
|
|
let octant_index = Octree::get_id_from_center(center_point, particle.position);
|
|
|
|
particle_octs[octant_index].push(particle);
|
|
|
|
@ -129,7 +139,7 @@ impl Octree { |
|
|
|
//Recurse
|
|
|
|
for (i, particle_oct) in particle_octs.iter().enumerate() {
|
|
|
|
self.add_particles_to_specific_node(particle_oct.clone(),
|
|
|
|
self.nodes[cur_node].children[i] as usize);
|
|
|
|
self.children[cur_node][i] as usize);
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
@ -179,18 +189,3 @@ impl Octree { |
|
|
|
(min_coord, max_coord)
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
impl OctNode {
|
|
|
|
fn new(min: Vector3D, max: Vector3D, parent: isize) -> Self {
|
|
|
|
OctNode {
|
|
|
|
parent: parent,
|
|
|
|
children: [-1, -1, -1, -1, -1, -1, -1, -1],
|
|
|
|
mass: 0.0,
|
|
|
|
min_coord: min,
|
|
|
|
max_coord: max,
|
|
|
|
num_particles: 0,
|
|
|
|
avg_weighted_position: Vector3D::new(0.0, 0.0, 0.0),
|
|
|
|
center_of_mass: Vector3D::new(0.0, 0.0, 0.0)
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
