My Sketch_tsa

/**

 * Java implementation of the Navier-Stokes-Solver from

 * http://www.dgp.toronto.edu/people/stam/reality/Research/pdf/GDC03.pdf

這個程序使用Navier-Stokes方程的離散化版本，通過求解一系列更新速度和密度的步驟，來模擬流體的行為。

1.密度和速度的初始化: 初始化流體密度場和速度場，以及相關的參數和初始條件。

2.Apply Forces (施加力): 根據滑鼠位置的變化施加外力（用戶互動），影響流體速度場。

3.Advect (移動): 將密度場和速度場根據速度場的流動進行移動。

4.Diffuse (擴散): 擴散速度場和密度場，模擬流體的黏性和擴散特性。

5.Project (投影): 計算速度場的投影，以保證速度場符合質量守恆和無旋流的特性。

 */

import processing.core.PVector;

​

float angleX = 0;

float angleY = 0;

float radius = 200;

​

NavierStokesSolver fluidSolver;

double visc, diff, vScale, velocityScale;

float  limitVelocity;

float oldMouseX = 1, oldMouseY = 1;

int numParticles;

Particle[] particles;

int c = color(255);

boolean vectors = false;

float zoom = 1.0;

​

​

public void setup() {

  //size(512, 512, P3D);

  size(512, 512);

  background(80, 0, 0);

  fill(150,0,0);

  fluidSolver = new NavierStokesSolver();

  numParticles = (int)pow(2, 16);

  particles = new Particle[numParticles];

  visc = 0.00025f;

  diff = 0.03f;

  velocityScale = 10;

  vectors = true;

  limitVelocity =10;

  initParticles();

}

​

private void initParticles() {

  for (int i = 0; i < numParticles - 1; i++) {

    particles[i] = new Particle();

    particles[i].x = random(width);

    particles[i].y = random(height);

  }

}

​

public void draw() {

  background(80, 0, 0);

  Camera();

  handleMouseMotion();

  //rect(0, 0, width, height);

  double dt = 1 / frameRate;

  fluidSolver.tick(dt, visc, diff);

  vScale = velocityScale * 60. / frameRate;

  drawParticlesPixels();

  // Draw particles in 3D space

  //drawParticles3D();

}

void Camera(){

 lights();  // Enable lighting

 if (mousePressed && key == CODED && keyCode == CONTROL) { //ctrl+left

    float camX = radius * cos(angleY) * sin(angleX);

    float camY = radius * sin(angleY);

    float camZ = radius * cos(angleY) * cos(angleX);

​

    camera(camX, camY, camZ, 0, 0, 0, 0, 1, 0);

    angleX += 0.01;

    angleY += 0.02;

  }

}

 void mouseWheel(MouseEvent event) {

  // 根据滚轮事件调整缩放

  //p5.js

  //zoom += event.delta / 500;

  zoom += event.getCount() / 500.0;

  // 限制缩放范围

  zoom = constrain(zoom, 0.1, 5);

}

​

private void drawParticlesPixels() {

  int n = NavierStokesSolver.N;

  float cellHeight = height / n;

  float cellWidth = width / n;

  for (Particle p : particles) {

    if (p != null) {

      int cellX = floor(p.x / cellWidth);

      int cellY = floor(p.y / cellHeight);

      float dx = (float) fluidSolver.getDx(cellX, cellY);

      float dy = (float) fluidSolver.getDy(cellX, cellY);

      float lX = p.x - cellX * cellWidth - cellWidth / 2;

      float lY = p.y - cellY * cellHeight - cellHeight / 2;

      int v, h, vf, hf;

      if (lX > 0) {

        v = Math.min(n, cellX + 1);

        vf = 1;

      } 

      else {

        v = Math.max(0, cellX - 1);

        vf = -1;

      }

      if (lY > 0) {

        h = Math.min(n, cellY + 1);

        hf = 1;

      } 

      else {

        h = Math.max(0, cellY - 1);

        hf = -1;

      }

      float dxv = (float) fluidSolver.getDx(v, cellY);

      float dxh = (float) fluidSolver.getDx(cellX, h);

      float dxvh = (float) fluidSolver.getDx(v, h);

      float dyv = (float) fluidSolver.getDy(v, cellY);

      float dyh = (float) fluidSolver.getDy(cellX, h);

      float dyvh = (float) fluidSolver.getDy(v, h);

      dx = lerp(lerp(dx, dxv, vf * lX / cellWidth), 

      lerp(dxh, dxvh, vf * lX / cellWidth), 

      hf * lY / cellHeight);

      dy = lerp(lerp(dy, dyv, vf * lX / cellWidth), 

      lerp(dyh, dyvh, vf * lX / cellWidth), 

      hf * lY / cellHeight);

      p.x += dx * vScale;

      p.y += dy * vScale;

      if (p.x < 0 || p.x >= width) {

        p.x = random(width);

      }

      if (p.y < 0 || p.y >= height) {

        p.y = random(height);

      }

      set((int) p.x, (int) p.y, c);

    }

  }

}

​

​

private void handleMouseMotion() {

  mouseX = max(1, mouseX);

  mouseY = max(1, mouseY);

  int n = NavierStokesSolver.N;

  float cellHeight = height / n;

  float cellWidth = width / n;

  float mouseDx = mouseX - oldMouseX;

  float mouseDy = mouseY - oldMouseY;

  int cellX = floor(mouseX / cellWidth);

  int cellY = floor(mouseY / cellHeight);

 // mouseDx = (abs(mouseDx) > limitVelocity) ? Math.signum(mouseDx) * limitVelocity : mouseDx;

 // mouseDy = (abs(mouseDy) > limitVelocity) ? Math.signum(mouseDy) * limitVelocity : mouseDy;

  fluidSolver.applyForce(cellX, cellY, mouseDx, mouseDy);

  oldMouseX = mouseX;

  oldMouseY = mouseY;

}

​

public class Particle {

  public float x;

  public float y;

​

}

​

​

public class NavierStokesSolver {

  final static int N = 16;

  final static int SIZE = (N + 2) * (N + 2);

  double[] u = new double[SIZE];

  double[] v = new double[SIZE];

  double[] u_prev = new double[SIZE];

  double[] v_prev = new double[SIZE];

  //  double[] dense = new double[SIZE];

  //  double[] dense_prev = new double[SIZE];

​

  public NavierStokesSolver() {

  }

​

  public double getDx(int x, int y) {

    return u[INDEX(x + 1, y + 1)];

  }

​

  public double getDy(int x, int y) {

    return v[INDEX(x + 1, y + 1)];

  }

​

  public void applyForce(int cellX, int cellY, double vx, double vy) {

    cellX += 1;

    cellY += 1;

    double dx = u[INDEX(cellX, cellY)];

    double dy = v[INDEX(cellX, cellY)];

​

    u[INDEX(cellX, cellY)] = (vx != 0) ? lerp((float) vx, 

      (float) dx, 0.85f) : dx;

    v[INDEX(cellX, cellY)] = (vy != 0) ? lerp((float) vy, 

      (float) dy, 0.85f) : dy;

  }

​

  void tick(double dt, double visc, double diff) {

    vel_step(u, v, u_prev, v_prev, visc, dt);

  }

​

  // method used to be 'static' since this class is not a top level type

  final int INDEX(int i, int j) {

    return i + (N + 2) * j;

  }

​

  double[] tmp = new double[SIZE];

​

  // same applies to the swap operation ^^

  final void SWAP(double[] x0, double[] x) {

   arraycopy(x0, 0, tmp, 0, SIZE);

   arraycopy(x, 0, x0, 0, SIZE);

   arraycopy(tmp, 0, x, 0, SIZE);

  }

​

  void add_source(double[] x, double[] s, double dt) {

    int i, size = (N + 2) * (N + 2);

    for (i = 0; i < size; i++)

      x[i] += dt * s[i];

  }

​

  void diffuse(int b, double[] x, double[] x0, double diff, double dt) {

    int i, j, k;

    double a = dt * diff * N * N;

    for (k = 0; k < 20; k++) {

      for (i = 1; i <= N; i++) {

        for (j = 1; j <= N; j++) {

          x[INDEX(i, j)] = (x0[INDEX(i, j)] + a

            * (x[INDEX(i - 1, j)] + x[INDEX(i + 1, j)]

            + x[INDEX(i, j - 1)] + x[INDEX(i, j + 1)]))

            / (1 + 4 * a);

        }

      }

      set_bnd(b, x);

    }

  }

​

  void advect(int b, double[] d, double[] d0, double[] u, double[] v, 

    double dt) {

    int i, j, i0, j0, i1, j1;

    double x, y, s0, t0, s1, t1, dt0;

    dt0 = dt * N;

    for (i = 1; i <= N; i++) {

      for (j = 1; j <= N; j++) {

        x = i - dt0 * u[INDEX(i, j)];

        y = j - dt0 * v[INDEX(i, j)];

        if (x < 0.5)

          x = 0.5;

        if (x > N + 0.5)

          x = N + 0.5;

        i0 = (int) x;

        i1 = i0 + 1;

        if (y < 0.5)

          y = 0.5;

        if (y > N + 0.5)

          y = N + 0.5;

        j0 = (int) y;

        j1 = j0 + 1;

        s1 = x - i0;

        s0 = 1 - s1;

        t1 = y - j0;

        t0 = 1 - t1;

        d[INDEX(i, j)] = s0

          * (t0 * d0[INDEX(i0, j0)] + t1 * d0[INDEX(i0, j1)])

          + s1

          * (t0 * d0[INDEX(i1, j0)] + t1 * d0[INDEX(i1, j1)]);

      }

    }

    set_bnd(b, d);

  }

​

  void set_bnd(int b, double[] x) {

    int i;

    for (i = 1; i <= N; i++) {

      x[INDEX(0, i)] = (b == 1) ? -x[INDEX(1, i)] : x[INDEX(1, i)];

      x[INDEX(N + 1, i)] = b == 1 ? -x[INDEX(N, i)] : x[INDEX(N, i)];

      x[INDEX(i, 0)] = b == 2 ? -x[INDEX(i, 1)] : x[INDEX(i, 1)];

      x[INDEX(i, N + 1)] = b == 2 ? -x[INDEX(i, N)] : x[INDEX(i, N)];

    }

    x[INDEX(0, 0)] = 0.5 * (x[INDEX(1, 0)] + x[INDEX(0, 1)]);

    x[INDEX(0, N + 1)] = 0.5 * (x[INDEX(1, N + 1)] + x[INDEX(0, N)]);

    x[INDEX(N + 1, 0)] = 0.5 * (x[INDEX(N, 0)] + x[INDEX(N + 1, 1)]);

    x[INDEX(N + 1, N + 1)] = 0.5 * (x[INDEX(N, N + 1)] + x[INDEX(N + 1, N)]);

  }

​

  void dens_step(double[] x, double[] x0, double[] u, double[] v, 

    double diff, double dt) {

    add_source(x, x0, dt);

    SWAP(x0, x);

    diffuse(0, x, x0, diff, dt);

    SWAP(x0, x);

    advect(0, x, x0, u, v, dt);

  }

​

  void vel_step(double[] u, double[] v, double[] u0, double[] v0, 

    double visc, double dt) {

    diffuse(1, u, u, visc, dt);

    diffuse(2, v, v, visc, dt);

    project(u, v, u0, v0);

  }

​

  void project(double[] u, double[] v, double[] p, double[] div) {

    int i, j, k;

    double h;

    h = 1.0 / N;

    for (i = 1; i <= N; i++) {

      for (j = 1; j <= N; j++) {

        div[INDEX(i, j)] = -0.5

          * h

          * (u[INDEX(i + 1, j)] - u[INDEX(i - 1, j)]

          + v[INDEX(i, j + 1)] - v[INDEX(i, j - 1)]);

        p[INDEX(i, j)] = 0;

      }

    }

    set_bnd(0, div);

    set_bnd(0, p);

    for (k = 0; k < 20; k++) {

      for (i = 1; i <= N; i++) {

        for (j = 1; j <= N; j++) {

          p[INDEX(i, j)] = (div[INDEX(i, j)] + p[INDEX(i - 1, j)]

            + p[INDEX(i + 1, j)] + p[INDEX(i, j - 1)] + p[INDEX(

            i, j + 1)]) / 4;

        }

      }

      set_bnd(0, p);

    }

    for (i = 1; i <= N; i++) {

      for (j = 1; j <= N; j++) {

        u[INDEX(i, j)] -= 0.5

          * (p[INDEX(i + 1, j)] - p[INDEX(i - 1, j)]) / h;

        v[INDEX(i, j)] -= 0.5

          * (p[INDEX(i, j + 1)] - p[INDEX(i, j - 1)]) / h;

      }

    }

    set_bnd(1, u);

    set_bnd(2, v);

  }

}

​

​

​