diff --git a/src/kd_tree.rs b/src/kd_tree.rs index c2aa3bb..138d9d8 100644 --- a/src/kd_tree.rs +++ b/src/kd_tree.rs @@ -1,17 +1,52 @@ -use super::gen_grid2d; +#[derive(Debug, Clone, PartialEq, PartialOrd)] +pub struct Grid2D { + pub x: f64, + pub y: f64, +} + +#[derive(Debug, Clone)] +pub struct Points2D { + pub points: Vec, +} #[allow(dead_code)] #[derive(Debug, Clone)] pub struct KDTree { id: usize, - position: gen_grid2d::Grid2D, + position: Grid2D, left: Option>, right: Option>, } +impl Grid2D { + #[allow(dead_code)] + pub fn new(x_: f64, y_: f64) -> Self { + Grid2D { x: x_, y: y_ } + } + + #[allow(dead_code)] + pub fn distance_square(&self, point: &Grid2D) -> f64 { + let dx2 = (self.x - point.x) * (self.x - point.x); + let dy2 = (self.y - point.y) * (self.y - point.y); + dx2 + dy2 + } +} + +impl Points2D { + #[allow(dead_code)] + pub fn new() -> Self { + Points2D { points: vec![] } + } + + #[allow(dead_code)] + pub fn push(&mut self, x_r: f64, y_r: f64) { + self.points.push(Grid2D { x: x_r, y: y_r }); + } +} + impl KDTree { #[allow(dead_code)] - fn new(vector: &gen_grid2d::Grid2D, id_: usize) -> Self { + fn new(vector: &Grid2D, id_: usize) -> Self { Self { id: id_, position: vector.clone(), @@ -43,20 +78,14 @@ impl KDTree { } #[allow(dead_code)] - pub fn neighbor_search(&self, x: &gen_grid2d::Grid2D, radius: f64) -> Vec { + pub fn neighbor_search(&self, x: &Grid2D, radius: f64) -> Vec { let mut near = vec![0; 0]; self.search_points_id(x, radius, &mut near, 0); near.clone() } #[allow(dead_code)] - fn search_points_id( - &self, - x: &gen_grid2d::Grid2D, - radius: f64, - near: &mut Vec, - mut depth: i32, - ) { + fn search_points_id(&self, x: &Grid2D, radius: f64, near: &mut Vec, mut depth: i32) { let axis = depth % 2; let r_self = self.position.distance_square(x).sqrt(); if r_self < radius { @@ -140,7 +169,7 @@ impl KDTree { } #[allow(dead_code)] - fn insert(&mut self, point: &gen_grid2d::Grid2D, mut depth: i32, id: usize) { + fn insert(&mut self, point: &Grid2D, mut depth: i32, id: usize) { let axis = depth % 2; match axis { @@ -198,7 +227,7 @@ impl KDTree { } #[allow(dead_code)] - fn create_kd_tree(&mut self, vec: &gen_grid2d::Points2D) -> Self { + fn create_kd_tree(&mut self, vec: &Points2D) -> Self { let depth = 0; for i in 1..vec.points.len() { self.insert(&vec.points[i], depth, i); @@ -217,11 +246,8 @@ impl KDTree { } #[allow(dead_code)] - pub fn construct_kd_tree(vec: &gen_grid2d::Points2D) -> KDTree { - let mut tree = KDTree::new( - &gen_grid2d::Grid2D::new(vec.points[0].x, vec.points[0].y), - 0, - ); + pub fn construct_kd_tree(vec: &Points2D) -> KDTree { + let mut tree = KDTree::new(&Grid2D::new(vec.points[0].x, vec.points[0].y), 0); tree.create_kd_tree(vec) } @@ -247,7 +273,7 @@ mod tests { let mut rng: rand::rngs::StdRng = rand::SeedableRng::from_seed(seed); let num_point: usize = 600; - let mut vec = gen_grid2d::Points2D::new(); + let mut vec = Points2D::new(); for _ in 0..num_point { let x_r = 2.0 * (rng.gen::() - 0.5); @@ -268,7 +294,7 @@ mod tests { let mut rng: rand::rngs::StdRng = rand::SeedableRng::from_seed(seed); let num_point: usize = 10; - let mut vec = gen_grid2d::Points2D::new(); + let mut vec = Points2D::new(); for _ in 0..num_point { let x_r = 2.0 * (rng.gen::() - 0.5); @@ -277,7 +303,7 @@ mod tests { } let tree = KDTree::construct_kd_tree(&mut vec); - let center = gen_grid2d::Grid2D { x: 0.0, y: 0.0 }; + let center = Grid2D { x: 0.0, y: 0.0 }; let radius = 0.4; let near = tree.neighbor_search(¢er, radius); @@ -291,7 +317,7 @@ mod tests { let mut rng: rand::rngs::StdRng = rand::SeedableRng::from_seed(seed); let num_point: usize = 10; - let mut vec = gen_grid2d::Points2D::new(); + let mut vec = Points2D::new(); for _ in 0..num_point { let x_r = 2.0 * (rng.gen::() - 0.5); @@ -300,7 +326,7 @@ mod tests { } let tree = KDTree::construct_kd_tree(&mut vec); - let center = gen_grid2d::Grid2D { x: 0.4, y: 0.3 }; + let center = Grid2D { x: 0.4, y: 0.3 }; let radius = 0.5; let near = tree.neighbor_search(¢er, radius); @@ -317,7 +343,7 @@ mod tests { let mut rng: rand::rngs::StdRng = rand::SeedableRng::from_seed(seed); let num_point: usize = 10; - let mut vec = gen_grid2d::Points2D::new(); + let mut vec = Points2D::new(); for _ in 0..num_point { let x_r = 2.0 * (rng.gen::() - 0.5); @@ -327,7 +353,7 @@ mod tests { let tree = KDTree::construct_kd_tree(&mut vec); - let center = gen_grid2d::Grid2D { x: 0.4, y: 0.3 }; + let center = Grid2D { x: 0.4, y: 0.3 }; let radius = 0.5; let mut near = vec![0 as usize; 0]; tree.search_points_id(¢er, radius, &mut near, 0); @@ -342,7 +368,7 @@ mod tests { let mut rnd = rand::thread_rng(); for _ in 0..10 { let num_point: i32 = rnd.gen_range(0..100000); - let mut vec = gen_grid2d::Points2D::new(); + let mut vec = Points2D::new(); for _ in 0..num_point { let x_r = rnd.gen::() as f64; let y_r = rnd.gen::() as f64; @@ -355,13 +381,13 @@ mod tests { } #[allow(dead_code)] - fn benchmark_pre(size: i32) -> gen_grid2d::Points2D { + fn benchmark_pre(size: i32) -> Points2D { use rand::Rng; let mut rnd = rand::thread_rng(); let max = size; let num_point: i32 = rnd.gen_range(0..max); - let mut vec = gen_grid2d::Points2D::new(); + let mut vec = Points2D::new(); for _ in 0..num_point { let x_r = rnd.gen::() as f64 / max as f64; let y_r = rnd.gen::() as f64 / max as f64; @@ -371,9 +397,9 @@ mod tests { } #[allow(dead_code)] - fn benchmark(vec: &gen_grid2d::Points2D) { + fn benchmark(vec: &Points2D) { let tree = KDTree::construct_kd_tree(vec); - let center = gen_grid2d::Grid2D::new(0.0, 0.0); + let center = Grid2D::new(0.0, 0.0); tree.neighbor_search(¢er, 0.01); } diff --git a/src/gen_grid2d.rs b/src/lennard_jones_potential.rs similarity index 60% rename from src/gen_grid2d.rs rename to src/lennard_jones_potential.rs index 3c0af05..451426f 100644 --- a/src/gen_grid2d.rs +++ b/src/lennard_jones_potential.rs @@ -1,45 +1,10 @@ use super::kd_tree; -#[derive(Debug, Clone, PartialEq, PartialOrd)] -pub struct Grid2D { - pub x: f64, - pub y: f64, -} - -impl Grid2D { - #[allow(dead_code)] - pub fn new(x_: f64, y_: f64) -> Self { - Grid2D { x: x_, y: y_ } - } - - #[allow(dead_code)] - pub fn distance_square(&self, point: &Grid2D) -> f64 { - let dx2 = (self.x - point.x) * (self.x - point.x); - let dy2 = (self.y - point.y) * (self.y - point.y); - dx2 + dy2 - } -} - -#[derive(Debug, Clone)] -pub struct Points2D { - pub points: Vec, -} - -impl Points2D { - #[allow(dead_code)] - pub fn new() -> Self { - Points2D { points: vec![] } - } - - #[allow(dead_code)] - pub fn push(&mut self, x_r: f64, y_r: f64) { - self.points.push(Grid2D { x: x_r, y: y_r }); - } - +impl kd_tree::Points2D { #[allow(dead_code)] pub fn euler_step_by_near_points( &mut self, - boundary: &Points2D, + boundary: &kd_tree::Points2D, tree: &kd_tree::KDTree, radius: f64, ) { @@ -72,7 +37,7 @@ impl Points2D { } #[allow(dead_code)] - pub fn euler_step(&mut self, boundary: &Points2D) { + pub fn euler_step(&mut self, boundary: &kd_tree::Points2D) { let dt = 1.0e-3; for i in 0..self.points.len() { let mut x = self.points[i].x - dt * self.lennard_jones_potential_deriv(i, boundary).x; @@ -99,10 +64,10 @@ impl Points2D { pub fn lennard_jones_potential_deriv_by_near_points( &mut self, index: usize, - boundary: &Points2D, + boundary: &kd_tree::Points2D, tree: &kd_tree::KDTree, radius: f64, - ) -> Grid2D { + ) -> kd_tree::Grid2D { let epsilon: f64 = 1.0; let sigma: f64 = 0.9 * (2.0_f64).powf(-1.0 / 6.0) / (self.points.len() as f64).sqrt(); let mut f_x = 0.0; @@ -123,16 +88,8 @@ impl Points2D { let r8 = r4 * r4; let r13 = r8 * r4 * r; let r7 = r2 * r4 * r; - f_x += 4.0 - * epsilon - * (12.0 * sigma12 / r13 - - 6.0 * sigma6 / r7) - * dx; - f_y += 4.0 - * epsilon - * (12.0 * sigma12 / r13 - - 6.0 * sigma6 / r7) - * dy; + f_x += 4.0 * epsilon * (12.0 * sigma12 / r13 - 6.0 * sigma6 / r7) * dx; + f_y += 4.0 * epsilon * (12.0 * sigma12 / r13 - 6.0 * sigma6 / r7) * dy; } } for i in 0..boundary.points.len() { @@ -144,54 +101,53 @@ impl Points2D { let r8 = r4 * r4; let r13 = r8 * r4 * r; let r7 = r2 * r4 * r; - f_x += 4.0 - * epsilon - * (12.0 * sigma12 / r13 - 6.0 * sigma6 / r7) - * dx; - f_y += 4.0 - * epsilon - * (12.0 * sigma12 / r13 - 6.0 * sigma6 / r7) - * dy; + f_x += 4.0 * epsilon * (12.0 * sigma12 / r13 - 6.0 * sigma6 / r7) * dx; + f_y += 4.0 * epsilon * (12.0 * sigma12 / r13 - 6.0 * sigma6 / r7) * dy; } - Grid2D { x: f_x, y: f_y } + kd_tree::Grid2D { x: f_x, y: f_y } } #[allow(dead_code)] - pub fn lennard_jones_potential_deriv(&mut self, index: usize, boundary: &Points2D) -> Grid2D { + pub fn lennard_jones_potential_deriv( + &mut self, + index: usize, + boundary: &kd_tree::Points2D, + ) -> kd_tree::Grid2D { let epsilon: f64 = 1.0; - let sigma: f64 = 0.9 * (2.0_f64).powf(-1.0 / 6.0) / (self.points.len() as f64).powf(0.5); + let sigma: f64 = 0.9 * (2.0_f64).powf(-1.0 / 6.0) / (self.points.len() as f64).sqrt(); let mut f_x = 0.0; let mut f_y = 0.0; + let sigma2 = sigma * sigma; + let sigma4 = sigma2 * sigma2; + let sigma8 = sigma4 * sigma4; + let sigma6 = sigma2 * sigma4; + let sigma12 = sigma8 * sigma4; for i in 0..self.points.len() { if i != index { let dx = self.points[index].x - self.points[i].x; let dy = self.points[index].y - self.points[i].y; - let r = (dx.powf(2.0) + dy.powf(2.0)).powf(0.5); - f_x += 4.0 - * epsilon - * (12.0 * sigma.powf(12.0) / r.powf(13.0) - - 6.0 * sigma.powf(6.0) / r.powf(7.0)) - * dx; - f_y += 4.0 - * epsilon - * (12.0 * sigma.powf(12.0) / r.powf(13.0) - - 6.0 * sigma.powf(6.0) / r.powf(7.0)) - * dy; + let r = (dx * dx + dy * dy).sqrt(); + let r2 = r * r; + let r4 = r2 * r2; + let r8 = r4 * r4; + let r13 = r8 * r4 * r; + let r7 = r2 * r4 * r; + f_x += 4.0 * epsilon * (12.0 * sigma12 / r13 - 6.0 * sigma6 / r7) * dx; + f_y += 4.0 * epsilon * (12.0 * sigma12 / r13 - 6.0 * sigma6 / r7) * dy; } } for i in 0..boundary.points.len() { let dx = self.points[index].x - boundary.points[i].x; let dy = self.points[index].y - boundary.points[i].y; - let r = (dx.powf(2.0) + dy.powf(2.0)).powf(0.5); - f_x += 4.0 - * epsilon - * (12.0 * sigma.powf(12.0) / r.powf(13.0) - 6.0 * sigma.powf(6.0) / r.powf(7.0)) - * dx; - f_y += 4.0 - * epsilon - * (12.0 * sigma.powf(12.0) / r.powf(13.0) - 6.0 * sigma.powf(6.0) / r.powf(7.0)) - * dy; + let r = (dx * dx + dy * dy).sqrt(); + let r2 = r * r; + let r4 = r2 * r2; + let r8 = r4 * r4; + let r13 = r8 * r4 * r; + let r7 = r2 * r4 * r; + f_x += 4.0 * epsilon * (12.0 * sigma12 / r13 - 6.0 * sigma6 / r7) * dx; + f_y += 4.0 * epsilon * (12.0 * sigma12 / r13 - 6.0 * sigma6 / r7) * dy; } - Grid2D { x: f_x, y: f_y } + kd_tree::Grid2D { x: f_x, y: f_y } } } diff --git a/src/main.rs b/src/main.rs index c1e4094..7566d39 100644 --- a/src/main.rs +++ b/src/main.rs @@ -5,11 +5,11 @@ use std::fs::File; use std::io::{BufWriter, Read}; use std::path::Path; -mod gen_grid2d; mod kd_tree; +mod lennard_jones_potential; #[allow(dead_code)] -fn draw_graph(i: usize, vec: &gen_grid2d::Points2D, boundary: &gen_grid2d::Points2D) { +fn draw_graph(i: usize, vec: &kd_tree::Points2D, boundary: &kd_tree::Points2D) { let out_file_name = format!("{:04}", i).to_string() + ".png"; let root = BitMapBackend::new(&out_file_name, (2000, 2000)).into_drawing_area(); @@ -93,8 +93,8 @@ fn main() { let num_point: usize = 100000; let num_boundary: usize = (std::f64::consts::PI * (num_point as f64).sqrt()) as usize; - let mut vec = gen_grid2d::Points2D::new(); - let mut boundary = gen_grid2d::Points2D::new(); + let mut vec = kd_tree::Points2D::new(); + let mut boundary = kd_tree::Points2D::new(); loop { let x_r = 2.0 * (rng.gen::() - 0.5); @@ -110,7 +110,7 @@ fn main() { let radius = 6.0 * std::f64::consts::PI / (num_boundary as f64); println!("radius = {}", radius); - let max_counter = 500; + let max_counter = 5_000; let mut tree = kd_tree::KDTree::construct_kd_tree(&vec);