inital commit

.gitignore

@@ -0,0 +1,2 @@
+/target
+**/*.rs.bk

Cargo.toml

@@ -0,0 +1,10 @@
+[package]
+name = "fluid"
+version = "0.1.0"
+authors = ["Cpt. Captain <nilskoch94@gmail.com>"]
+edition = "2018"
+
+[dependencies]
+piston = "0.42.0"
+piston_window = "0.89.0"
+piston2d-gfx_graphics = "0.56.0"

src/fluid.rs

@@ -0,0 +1,294 @@
+pub struct FluidSquare {
+    size: i32,
+    dt: f32,
+    diff: f32,
+    visc: f32,
+
+    s: Vec<f32>,
+    pub density: Vec<f32>,
+
+    pub v_x: Vec<f32>,
+    pub v_y: Vec<f32>,
+
+    v_x0: Vec<f32>,
+    v_y0: Vec<f32>,
+}
+
+impl FluidSquare {
+    pub    fn new(size: i32, diff: f32, visc: f32, dt: f32) -> FluidSquare {
+        FluidSquare {
+            size: size,
+            dt: dt,
+            diff: diff,
+            visc: visc,
+            s: vec![0.0; (size * size) as usize],
+            density: vec![0.0; (size * size) as usize],
+            v_x: vec![0.0; (size * size) as usize],
+            v_y: vec![0.0; (size * size) as usize],
+            v_x0: vec![0.0; (size * size) as usize],
+            v_y0: vec![0.0; (size * size) as usize],
+        }
+    }
+
+    pub fn step(&mut self, iter: i32) {
+        diffuse(
+            1,
+            &mut self.v_x0,
+            &mut self.v_x,
+            self.visc,
+            self.dt,
+            iter,
+            self.size,
+        );
+        diffuse(
+            2,
+            &mut self.v_y0,
+            &mut self.v_y,
+            self.visc,
+            self.dt,
+            iter,
+            self.size,
+        );
+
+        project(
+            &mut self.v_x0,
+            &mut self.v_y0,
+            &mut self.v_x,
+            &mut self.v_y,
+            iter,
+            self.size,
+        );
+
+        advect(
+            1,
+            &mut self.v_x,
+            &self.v_x0,
+            &self.v_x0,
+            &self.v_y0,
+            self.dt,
+            self.size,
+        );
+        advect(
+            2,
+            &mut self.v_y,
+            &self.v_y0,
+            &self.v_x0,
+            &self.v_y0,
+            self.dt,
+            self.size,
+        );
+
+        project(
+            &mut self.v_x,
+            &mut self.v_y,
+            &mut self.v_x0,
+            &mut self.v_y0,
+            iter,
+            self.size,
+        );
+
+        diffuse(
+            0,
+            &mut self.s,
+            &mut self.density,
+            self.diff,
+            self.dt,
+            iter,
+            self.size,
+        );
+        advect(
+            0,
+            &mut self.density,
+            &self.s,
+            &self.v_x,
+            &self.v_y,
+            self.dt,
+            self.size,
+        );
+    }
+
+    pub fn add_dye(&mut self, x: i32, y: i32, amount: f32) {
+        let index = ix(x, y, self.size);
+        self.density[index as usize] += amount;
+    }
+
+    pub fn add_velocity(&mut self, x: i32, y: i32, amount_x: f32, amount_y: f32) {
+        let index = ix(x, y, self.size);
+        self.v_x[index as usize] += amount_x;
+        self.v_y[index as usize] += amount_y;
+    }
+}
+
+fn diffuse(b: i32, x: &mut Vec<f32>, x0: &Vec<f32>, diff: f32, dt: f32, iter: i32, size: i32) {
+    let a: f32 = dt * diff * (size - 2) as f32 * (size - 2) as f32;
+    lin_solve(b, x, x0, a, 1.0 + 6.0 * a, iter, size)
+}
+
+fn lin_solve(b: i32, x: &mut Vec<f32>, x0: &Vec<f32>, a: f32, c: f32, iter: i32, size: i32) {
+    let c_recip = 1.0 / c;
+    for _k in 0..iter {
+        for j in 1..size - 1 {
+            for i in 1..size - 1 {
+                x[ix(i, j, size) as usize] = (x0[ix(i, j, size) as usize]
+                    + a * (x[ix(i + 1, j, size) as usize]
+                        + x[ix(i - 1, j, size) as usize]
+                        + x[ix(i, j + 1, size) as usize]
+                        + x[ix(i, j - 1, size) as usize]))
+                    * c_recip;
+            }
+        }
+    }
+    set_bnd(b, x, size);
+}
+
+fn set_bnd(b: i32, x: &mut Vec<f32>, size: i32) {
+    for i in 1..size - 1 {
+        x[ix(i, 0, size) as usize] = if b == 2 {
+            -x[ix(i, 1, size) as usize]
+        } else {
+            x[ix(i, 1, size) as usize]
+        };
+        x[ix(i, size - 1, size) as usize] = if b == 2 {
+            -x[ix(i, size - 2, size) as usize]
+        } else {
+            x[ix(i, size - 2, size) as usize]
+        };
+    }
+
+    for j in 1..size - 1 {
+        x[ix(0, j, size) as usize] = if b == 1 {
+            -x[ix(1, j, size) as usize]
+        } else {
+            x[ix(1, j, size) as usize]
+        };
+        x[ix(size - 1, j, size) as usize] = if b == 1 {
+            -x[ix(size - 2, j, size) as usize]
+        } else {
+            x[ix(size - 2, j, size) as usize]
+        };
+    }
+
+    x[ix(0, 0, size) as usize] = 0.5 * (x[ix(1, 0, size) as usize] + x[ix(0, 1, size) as usize]);
+
+    x[ix(0, size - 1, size) as usize] =
+        0.5 * (x[ix(1, size - 1, size) as usize] + x[ix(0, size - 2, size) as usize]);
+    x[ix(size - 1, 0, size) as usize] =
+        0.5 * (x[ix(size - 2, 0, size) as usize] + x[ix(size - 1, 1, size) as usize]);
+    x[ix(size - 1, size - 1, size) as usize] =
+        0.5 * (x[ix(size - 2, size - 1, size) as usize] + x[ix(size - 1, size - 2, size) as usize]);
+}
+
+fn project(
+    v_x: &mut Vec<f32>,
+    v_y: &mut Vec<f32>,
+    p: &mut Vec<f32>,
+    div: &mut Vec<f32>,
+    iter: i32,
+    size: i32,
+) {
+    for j in 1..size - 1 {
+        for i in 1..size - 1 {
+            div[ix(i, j, size) as usize] = -0.5
+                * (v_x[ix(i + 1, j, size) as usize] - v_x[ix(i - 1, j, size) as usize]
+                    + v_y[ix(i, j + 1, size) as usize]
+                    - v_y[ix(i, j - 1, size) as usize])
+                / size as f32;
+            p[ix(i, j, size) as usize] = 0.0;
+        }
+    }
+    set_bnd(0, div, size);
+    set_bnd(0, p, size);
+    lin_solve(0, p, div, 1.0, 6.0, iter, size);
+
+    for j in 1..size - 1 {
+        for i in 1..size - 1 {
+            v_x[ix(i, j, size) as usize] -= 0.5
+                * (p[ix(i + 1, j, size) as usize] - p[ix(i - 1, j, size) as usize])
+                * size as f32;
+            v_y[ix(i, j, size) as usize] -= 0.5
+                * (p[ix(i, j + 1, size) as usize] - p[ix(i, j - 1, size) as usize])
+                * size as f32;
+        }
+    }
+    set_bnd(1, v_x, size);
+    set_bnd(2, v_y, size);
+}
+
+fn advect(
+    b: i32,
+    d: &mut Vec<f32>,
+    d0: &Vec<f32>,
+    v_x: &Vec<f32>,
+    v_y: &Vec<f32>,
+    dt: f32,
+    size: i32,
+) {
+    let (mut i0, mut i1, mut j0, mut j1): (f32, f32, f32, f32);
+    let dtx: f32 = dt * (size - 2) as f32;
+    let dty: f32 = dt * (size - 2) as f32;
+
+    let (mut s0, mut s1, mut t0, mut t1): (f32, f32, f32, f32);
+    let (mut tmp1, mut tmp2, mut x, mut y): (f32, f32, f32, f32);
+
+    let nfloat: f32 = size as f32;
+    let (mut ifloat, mut jfloat): (f32, f32);
+
+    for i in 1..size - 1 {
+        ifloat = i as f32;
+        for j in 1..size - 1 {
+            jfloat = j as f32;
+
+            tmp1 = dtx * v_x[ix(i, j, size) as usize];
+            tmp2 = dty * v_y[ix(i, j, size) as usize];
+            x = ifloat - tmp1;
+            y = jfloat - tmp2;
+
+            if x < 0.5 {
+                x = 0.5
+            }
+            if x > nfloat + 0.5 {
+                x = nfloat + 0.5
+            }
+            i0 = x.floor();
+            i1 = i0 + 1.0;
+
+            if y < 0.5 {
+                y = 0.5
+            }
+            if y > nfloat + 0.5 {
+                y = nfloat + 0.5
+            }
+            j0 = y.floor();
+            j1 = j0 + 1.0;
+
+            s1 = x - i0;
+            s0 = 1.0 - s1;
+            t1 = y - j0;
+            t0 = 1.0 - t1;
+
+            let i0i: i32 = i0 as i32;
+            let i1i: i32 = i1 as i32;
+            let j0i: i32 = j0 as i32;
+            let j1i: i32 = j1 as i32;
+
+            // probably wrong
+                        d[ix(i, j, size) as usize] = s0
+                            * (t0 * d0[ix(i0i, j0i, size) as usize] + t1 * d0[ix(i0i, j1i, size) as usize])
+                            + s1 * (t0 * d0[ix(i1i, j0i, size) as usize] + t1 * d0[ix(i1i, j1i, size) as usize])
+
+//             d[ix(i, j, size) as usize] = s0
+//                 * (t0 * d0[ix(i0i, j0i, size) as usize] + t1 * d0[ix(i0i, j1i, size) as usize])
+//                 + s1 * (t0 * d0[ix(i1i, j0i, size) as usize]
+//                     + t1 * d0[ix(i1i, j1i, size) as usize]);
+        }
+    }
+
+    set_bnd(b, d, size);
+}
+
+pub fn ix(x: i32, y: i32, size: i32) -> i32 {
+    // (x + (y * size)).max(0).min(size * size - 1)
+    let r_x = x.max(0).min(size - 1);
+    let r_y = y.max(0).min(size - 1);
+    r_x + (r_y * size)
+}

src/main.rs

@@ -0,0 +1,120 @@
+/* A Rust implementation of Mike Ash's fluid solver from
+ * https://mikeash.com/pyblog/fluid-simulation-for-dummies.html
+ *
+ * Inspired by The Coding Trains Challenge #132
+ * https://www.youtube.com/watch?v=alhpH6ECFvQ
+ */
+
+extern crate piston_window;
+
+use piston_window::*;
+
+mod fluid;
+mod tests;
+
+fn main() {
+    let size: i32 = 80;
+    let scale: i32 = 5;
+    let iter: i32 = 10;
+    let dye_amount = 1.0;
+    let diffusion = 0.2;
+    let viscosity = 0.000000001;
+    let dt = 0.000001;
+    let mut fs: fluid::FluidSquare = fluid::FluidSquare::new(size, dt, viscosity, diffusion);
+    // fs.step(iter);
+    println!("Hello, world!");
+
+    let _background = [0.1, 0.1, 0.1, 1.0];
+    let mut events = Events::new(EventSettings::new());
+    let mut dye = false;
+    let mut vel = false;
+    let mut vis = true;
+    let mut prev_x: f64 = 0.0;
+    let mut prev_y: f64 = 0.0;
+
+    let mut window: PistonWindow = WindowSettings::new("Hello Piston!", [(size * scale) as u32; 2])
+        .exit_on_esc(true)
+        // .fullscreen(true)
+        .build()
+        .unwrap_or_else(|e| panic!("Failed to build PistonWindow: {}", e));
+
+    while let Some(e) = events.next(&mut window) {
+        window.draw_2d(&e, |_c, g| {
+            fs.step(iter);
+            clear(_background, g);
+            //  clear(color::BLACK, g);
+            for i in 0..size {
+                for j in 0..size {
+                    let x = i as i32 * scale;
+                    let y = j as i32 * scale;
+
+                    let rectangle = Rectangle::new(color::alpha(if vis {
+                        fs.density[fluid::ix(i, j, size) as usize] // - 0.02,
+                    } else {
+                        0.5 * (fs.v_x[fluid::ix(i, j, size) as usize]
+                            + fs.v_y[fluid::ix(i, j, size) as usize])
+                            .abs()
+                    }));
+
+                    // let square = rectangle::square(
+                    //                         x as f64,
+                    //                         y as f64,
+                    //                         scale as f64,
+                    //                         )
+                    let square =
+                        rectangle::centered([x as f64, y as f64, scale as f64, scale as f64]);
+                    let transform = _c.transform.trans(0.0, 0.0);
+
+                    rectangle.draw(square, &_c.draw_state, transform, g);
+                }
+            }
+        });
+        if let Some(p) = e.press_args() {
+            if p == Button::Mouse(MouseButton::Left) {
+                dye = true;
+                vel = true;
+            }
+            if p == Button::Mouse(MouseButton::Right) {
+                vel = true;
+            }
+            if p == Button::Keyboard(Key::V) {
+                vis = !vis;
+            }
+        }
+        if let Some(p) = e.release_args() {
+            if p == Button::Mouse(MouseButton::Left) {
+                dye = false;
+                vel = false;
+            }
+            if p == Button::Mouse(MouseButton::Right) {
+                vel = false;
+            }
+        }
+
+        if let Some(m) = e.mouse_cursor_args() {
+            if dye {
+                fs.add_dye(
+                    //                  (e[0] as i32 / (scale * scale)) as i32,
+                    //                  (e[1] as i32 / (scale * scale)) as i32,
+                    (m[0] as i32 / scale) as i32,
+                    (m[1] as i32 / scale) as i32,
+                    dye_amount,
+                );
+            }
+            if vel {
+                let mut vel_add = 50.0;
+                if dye {
+                    vel_add = 20.0;
+                }
+                fs.add_velocity(
+                    (m[0] as i32 / scale) as i32,
+                    (m[1] as i32 / scale) as i32,
+                    ((m[0] - prev_x) * vel_add) as f32,
+                    ((m[1] - prev_y) * vel_add) as f32,
+                );
+                prev_x = m[0];
+                prev_y = m[1];
+            }
+        }
+    }
+}