Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

高速化 #5

Merged
merged 1 commit into from
Aug 8, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 3 additions & 0 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -16,4 +16,7 @@ Cargo.lock
*.png
*.log

*.svg
perf.*

!images/*.png
37 changes: 26 additions & 11 deletions src/gen_grid2d.rs
Original file line number Diff line number Diff line change
Expand Up @@ -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
}
}
Expand Down Expand Up @@ -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 }
Expand Down
16 changes: 8 additions & 8 deletions src/main.rs
Original file line number Diff line number Diff line change
Expand Up @@ -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();

Expand All @@ -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),
))
Expand All @@ -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),
))
Expand Down Expand Up @@ -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();
Expand All @@ -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);

Expand All @@ -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);
Expand Down