diff --git a/.gitignore b/.gitignore index 54170cf..c9c4707 100644 --- a/.gitignore +++ b/.gitignore @@ -16,4 +16,7 @@ Cargo.lock *.png *.log +*.svg +perf.* + !images/*.png \ No newline at end of file diff --git a/src/gen_grid2d.rs b/src/gen_grid2d.rs index 40051c0..3c0af05 100644 --- a/src/gen_grid2d.rs +++ b/src/gen_grid2d.rs @@ -14,8 +14,8 @@ impl Grid2D { #[allow(dead_code)] pub fn distance_square(&self, point: &Grid2D) -> f64 { - let dx2 = (self.x - point.x).powf(2.0); - let dy2 = (self.y - point.y).powf(2.0); + let dx2 = (self.x - point.x) * (self.x - point.x); + let dy2 = (self.y - point.y) * (self.y - point.y); dx2 + dy2 } } @@ -104,38 +104,53 @@ impl Points2D { radius: f64, ) -> 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; let near = tree.neighbor_search(&self.points[index], radius); for k in near.iter() { if index != *k { let dx = self.points[index].x - self.points[*k].x; let dy = self.points[index].y - self.points[*k].y; - let r = (dx.powf(2.0) + dy.powf(2.0)).powf(0.5); + 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 * sigma.powf(12.0) / r.powf(13.0) - - 6.0 * sigma.powf(6.0) / r.powf(7.0)) + * (12.0 * sigma12 / r13 + - 6.0 * sigma6 / r7) * 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)) + * (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); + 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 * sigma.powf(12.0) / r.powf(13.0) - 6.0 * sigma.powf(6.0) / r.powf(7.0)) + * (12.0 * sigma12 / r13 - 6.0 * sigma6 / r7) * 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)) + * (12.0 * sigma12 / r13 - 6.0 * sigma6 / r7) * dy; } Grid2D { x: f_x, y: f_y } diff --git a/src/main.rs b/src/main.rs index 69e7e0c..c1e4094 100644 --- a/src/main.rs +++ b/src/main.rs @@ -12,7 +12,7 @@ mod kd_tree; fn draw_graph(i: usize, vec: &gen_grid2d::Points2D, boundary: &gen_grid2d::Points2D) { let out_file_name = format!("{:04}", i).to_string() + ".png"; - let root = BitMapBackend::new(&out_file_name, (1080, 1080)).into_drawing_area(); + let root = BitMapBackend::new(&out_file_name, (2000, 2000)).into_drawing_area(); root.fill(&WHITE).unwrap(); @@ -29,7 +29,7 @@ fn draw_graph(i: usize, vec: &gen_grid2d::Points2D, boundary: &gen_grid2d::Point chart .draw_series(PointSeries::of_element( (0..vec.points.len()).map(|i| (vec.points[i].x, vec.points[i].y)), - 2, + 1, ShapeStyle::from(&RED).filled(), &|coord, size, style| EmptyElement::at(coord) + Circle::new((0, 0), size, style), )) @@ -38,7 +38,7 @@ fn draw_graph(i: usize, vec: &gen_grid2d::Points2D, boundary: &gen_grid2d::Point chart .draw_series(PointSeries::of_element( (0..boundary.points.len()).map(|i| (boundary.points[i].x, boundary.points[i].y)), - 2, + 1, ShapeStyle::from(&BLUE).filled(), &|coord, size, style| EmptyElement::at(coord) + Circle::new((0, 0), size, style), )) @@ -90,8 +90,8 @@ fn main() { let seed: [u8; 32] = [1; 32]; let mut rng: rand::rngs::StdRng = rand::SeedableRng::from_seed(seed); - let num_point: usize = 400; - let num_boundary: usize = (std::f64::consts::PI * (num_point as f64).powf(0.5)) as usize; + 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(); @@ -107,10 +107,10 @@ fn main() { } } - let radius = 4.0 * std::f64::consts::PI / (num_boundary as f64); + let radius = 6.0 * std::f64::consts::PI / (num_boundary as f64); println!("radius = {}", radius); - let max_counter = 50; + let max_counter = 500; let mut tree = kd_tree::KDTree::construct_kd_tree(&vec); @@ -121,7 +121,7 @@ fn main() { for i in 0..max_counter { draw_graph(i, &vec, &boundary); - for _ in 0..100 { + for _ in 0..10 { vec.euler_step_by_near_points(&boundary, &tree, radius); //vec.euler_step(&boundary); tree = kd_tree::KDTree::construct_kd_tree(&vec);