3D fluid simulation

prosession

Hey guys, I'm currently working on a project for simulating rigid bodies as fluids. The fluid acts inside the boundaries of an imported OBJ, like an aquarium. I found some nice papers about simulating fluids by Navier-Stokes Equations here and here. I also found a processing code for applying these equations in 2D but I still can't figure out how to add the third dimension. Honestly, I think the hole concept of the Navier-Stokes Equations isn't 100% clear to me. I'm grateful for any information about these algorithms that fill my gaps.

int n = 128;          // The size of the calculation grid
int gridSize = n+2;   // Extra grid space for boundary
int pixelSize = 2;    // The size of each grid square on the screen

// I unravelled the 1D arrays from Jos Stam's paper back to 2D arrays, as we don't have compile time macros in Java...
float[][] u = new float[gridSize][gridSize];
float[][] v = new float[gridSize][gridSize];
float[][] uPrev = new float[gridSize][gridSize];
float[][] vPrev = new float[gridSize][gridSize];
float[][] dens = new float[gridSize][gridSize];
float[][] densPrev = new float[gridSize][gridSize];

float viscosity = 0.0001;  // Viscosity of fluid
float dt = 0.2;            // Rate of change
float diff = 0.0001;       // Degree of diffusion of density over time

int prevMouseX;
int prevMouseY;

// Flags for control of display
boolean showDensity = true;
boolean showVelocity = false;
boolean showNormals = false;
boolean showLighting = true;
boolean showGradient = false;

// Flag for drawing velocity or density - shift (or RMB) = density
boolean shiftPressed = false;

Gradient gradient = new Gradient();
color[] gradientColours;
int numGradientColours = 256;  // Use a larger array for better gradient resolution

int densityBrushSize = n/10;   // Size of the density area applied with the mouse
int velocityBrushSize = n/20;  // Ditto velocity
int lineSpacing = n/20;        // Spacing between velocity and normal lines


void setup() {
  size ( n*pixelSize, n*pixelSize, OPENGL );
  background ( 0 );
  noStroke();

  gradient.makeDefaultGradient(gradient.DEFAULT_BLACKBODY);
  gradientColours = gradient.makeArrayOfColours(numGradientColours); 

  reset();

  frameRate(1000);
}

void draw() {
  setForce();
  calculateVelocity(u, v, uPrev, vPrev, viscosity, dt);
  calculateDensity(dens, densPrev, u, v, diff, dt);
  drawField(dens);
}

void reset() {
  initVelocity();
  initDensity();
}

void initField(float[][] f) {
  for ( int i=0; i < gridSize; i++ )
    for ( int j=0; j < gridSize; j++ )
      f[i][j] = 0.0;
}

void initVelocity() {
  initField(u);
  initField(v);
  initField(uPrev);
  initField(vPrev);
}

void initDensity() {
  initField(dens);
  initField(densPrev);
}

int ix(int i, int j) {
  return i + j*(n+2);
}

void addSource ( float[][] x, float[][] s, float dt )
{
  for ( int i=0; i<gridSize; i++ )
    for ( int j=0; j<gridSize; j++ )
      x[i][j] += s[i][j]*dt;
}

void setBnd ( int b, float[][] x )
{

  int i;
  for ( i=1; i <= n; i++ ) {

    if ( b==1 )  x [0  ][i  ] = -x [1  ][i  ];  else  x [0  ][i  ] = x [1  ][i  ];
    if ( b==1 )  x [n+1][i  ] = -x [n  ][i  ];  else  x [n+1][i  ] = x [n  ][i  ];
    if ( b==2 )  x [i  ][0  ] = -x [i  ][1  ];  else  x [i  ][0  ] = x [i  ][1  ];
    if ( b==2 )  x [i  ][n+1] = -x [i  ][n  ];  else  x [i  ][n+1] = x [i  ][n  ];

  }
  x [0  ][0  ] = 0.5 * ( x [1  ][0  ] + x [0  ][1  ] );
  x [0  ][n+1] = 0.5 * ( x [1  ][n+1] + x [0  ][n  ] );
  x [n+1][0  ] = 0.5 * ( x [n  ][0  ] + x [n+1][1  ] );
  x [n+1][n+1] = 0.5 * ( x [n  ][n+1] + x [n+1][n  ] );

}

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

  int i, j, k;
  float a = dt * diff * n * n;

  for ( k=0; k<20; k++ ) {
    for ( i=1; i <= n; i++ ) {
      for ( j=1; j <= n; j++ ) {
        x[i][j] = (x0[i][j] + a * ( x[i-1][j] + x[i+1][j] + x[i][j-1] + x[i][j+1] ) ) / (1+4*a);
      }
    }   
    setBnd(b, x);
  }

}

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

  int i, j, k;
  float h;

  h = 1.0/n;
  for ( i=1; i<=n; i++ ) {
     for ( j=1; j<=n; j++ ) {
        div [i][j] = -0.5 * h * ( u [i+1][j] - u [i-1][j] + v [i][j+1] - v [i][j-1] );
        p [i][j] = 0;
     }   
  }
  setBnd ( 0, div );
  setBnd ( 0, p );

  for ( k=0; k<20; k++ ) {
    for ( i=1; i<=n; i++ ) {
      for ( j=1; j<=n; j++ ) {
        p [i][j] = ( div [i][j] + p [i-1][j] + p [i+1][j] + p [i][j-1] + p [i][j+1] ) / 4;
      }
    }
    setBnd ( 0, p );
  }

  for ( i=1; i<=n; i++ ) {
    for ( j=1; j<=n; j++ ) {
      u [i][j] -= 0.5 * ( p [i+1][j] - p [i-1][j] ) / h;
      v [i][j] -= 0.5 * ( p [i][j+1] - p [i][j-1] ) / h;
    }   
  }
  setBnd ( 1, u );
  setBnd ( 2, v );

}

void advect ( int b, float[][] d, float[][] d0, float[][] u, float[][] v, float dt ) {

  int i, j, i0, j0, i1, j1;
  float 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 [i][j];
      y = j - dt0 * v [i][j];

      x = max (0.5, x);
      x = min (n+0.5, x);

      //i0 = (int)x;
      i0 = floor(x);
      i1 = i0 + 1;

      y = max (0.5, y);
      y = min (n+0.5, y);

      //j0 = (int)j;
      j0 = floor(y);
      j1 = j0+1;

      s1 = x-i0;
      s0 = 1-s1;
      t1 = y-j0;
      t0 = 1-t1;

      d [i][j] = s0 * ( t0 * d0 [i0][j0] + t1 * d0 [i0][j1] ) +
                    s1 * ( t0 * d0 [i1][j0] + t1 * d0 [i1][j1] );
    }
  }
  setBnd ( b, d );
}

void calculateVelocity(float[][] u, float[][] v, float[][] u0, float[][] v0, float visc, float dt) {

  addSource ( u, u0, dt );
  addSource ( v, v0, dt );

  float[][] tmp; 
  tmp = u;  u = u0;  u0 = tmp;
  tmp = v;  v = v0;  v0 = tmp;

  diffuse ( 1, u, u0, visc, dt );
  diffuse ( 2, v, v0, visc, dt );

  project ( u, v, u0, v0 );

  tmp = u;  u = u0;  u0 = tmp;
  tmp = v;  v = v0;  v0 = tmp;

  advect ( 1, u, u0, u0, v0, dt );
  advect ( 2, v, v0, u0, v0, dt );

  project ( u, v, u0, v0 );

}

void calculateDensity(float[][] x, float[][] x0, float[][] u, float[][] v, float diff, float dt) {
  float[][] tmp;

  addSource ( x, x0, dt );
  tmp = x; x = x0; x0 = tmp;
  diffuse ( 0, x, x0, diff, dt );
  tmp = x; x = x0; x0 = tmp;
  advect ( 0, x, x0, u, v, dt );
}

void drawField(float[][] f) {
  int x,y;
  float s;
  color col;

  float ax, ay, az;
  float bx, by, bz;
  float nx=0, ny=0, nz=0;

  float vu, vv;

  float l;

  float lx=0, ly=0, lz=0;
  float vx=0, vy=0, vz=0;
  float hx=0, hy=0, hz=0;

  float ndoth, ndotl;

  float w;
  float d;

  background(0);


  for (y = 1; y <= n; y++ ) {
    for ( x = 1; x <= n; x++ ) {

      d = dens[x][y];

      if ( showNormals || showLighting ) {
        az = d - dens[x+1][y];
        bz = d - dens[x][y+1];

        nx = -az;
        ny = -bz;
        nz = 1;

        l = -sqrt(nx*nx + ny*ny + 1); nx /= l; ny /= l; nz /= l;
      }

      if ( showDensity ) {
        s = abs(range(numGradientColours * d,0,numGradientColours-1));
        noStroke();

        col = gradientColours[int(s)];

        if ( showLighting ) {

          ndoth = nz;

          w = ((s*512)/numGradientColours) * pow(ndoth, 100000);

          col = blendColor ( col, color(w), ADD );

        }       

        fill ( col );
        rect ( (x-1) * pixelSize, (y-1) * pixelSize, pixelSize, pixelSize );
      }

      if ( showVelocity ) {
        if ( (x % lineSpacing) == 0 && (y % lineSpacing) == 0 ) {
          noFill();
          stroke(255,0,0);
          vu = range(500 * u[x][y], -50,50);
          vv = range(500 * v[x][y], -50,50);
          line( (x-1)*pixelSize, (y-1)*pixelSize, (x-1)*pixelSize + vu, (y-1)*pixelSize + vv);
        }
      }

      if ( showNormals ) {
        if ( (x % lineSpacing) == 0 && (y % lineSpacing) == 0 ) {
          noFill();
          stroke(255);       

          vu = range(500 * nx, -50,50);
          vv = range(500 * ny, -50,50);
          line( (x-1)*pixelSize, (y-1)*pixelSize, (x-1)*pixelSize + vu, (y-1)*pixelSize + vv);
        }
      }     
    }   
  }

  if ( showGradient ) {
    noStroke();
    for ( int i=0; i < numGradientColours; i++ ) {
      fill(gradientColours[i]);
      rect ( i % width, i / width, 1,10 );
    }
  }
}


void setForce() {
  initField(densPrev);
  initField(uPrev);
  initField(vPrev);

  if ( mousePressed ) {
    int x, y;

    x = ( mouseX * n ) / ( width ) + 1;
    y = ( mouseY * n ) / ( height ) + 1;

    if ( !keyPressed )
      shiftPressed = false;

    if ( mouseButton == LEFT && !shiftPressed ) {
      //uPrev [x,y] += (mouseX - prevMouseX);
      //vPrev [x,y] += (mouseY - prevMouseY);

      setForceArea ( uPrev, x, y, mouseX - prevMouseX, velocityBrushSize );
      setForceArea ( vPrev, x, y, mouseY - prevMouseY, velocityBrushSize );
    }
    else {
      //densPrev [x,y] += 5;
      int m = (mouseX - prevMouseX) + (mouseY - prevMouseY);
      setForceArea ( densPrev, x, y, range(abs(m),0,2), densityBrushSize );
    }
  }

  prevMouseX = mouseX;
  prevMouseY = mouseY;

}

float range(float f, float minf, float maxf) {
  return max ( min ( f, maxf ), minf );
}

void setForceArea(float[][] field, int x, int y, float s, float r) {

  int i,j, dx, dy;
  float f;

  for ( i = int(range(x-r,1,n)); i <= int(range(x+r,1,n)); i++ ) {
    dx = x - i;
    for ( j = int(range(y-r,1,n)); j <= int(range(y+r,1,n)); j++ ) {
      dy = y - j;
      f = 1 - ( sqrt(dx*dx + dy*dy) / r );
      field[i][j] += range(f,0,1) * s;
    }
  }

}


void keyPressed() {

  if ( keyCode == 32 ) {
    gradient.makeRandomGradient(4);
    gradientColours = gradient.makeArrayOfColours(numGradientColours);
  }

  if ( key == 'v' )
    showVelocity = !showVelocity;

  if ( key == 'n' )
    showNormals = !showNormals; 

  if ( key == 'd' )
    showDensity = !showDensity;

  if ( key == 'l' )
    showLighting = !showLighting;

  if ( key == 'r' )
    reset();

  if ( key == '+' )
    viscosity += 0.0001;
  if ( key == '-' )
    viscosity -= 0.0001;

  viscosity = max(viscosity, 0);

  if ( key >= '0' && key <= '7' ) {
    gradient.makeDefaultGradient(key - '0');
    gradientColours = gradient.makeArrayOfColours(numGradientColours);
  }

  if ( key == 'g' )
    showGradient = ! showGradient;

  if ( key == 's' )
    save("screenshot-" + year() + month() + day() + hour() + minute() + second() + ".png");

  if ( key == CODED )
    if ( keyCode == SHIFT )
      shiftPressed = true;

}

class Gradient {

  public static final int DEFAULT_BLACK_TO_WHITE = 0;
  public static final int DEFAULT_RANDOM = 1; // makes a new random gradient each time
  public static final int DEFAULT_SPECTRUM = 2;
  public static final int DEFAULT_INFRARED = 3;
  public static final int DEFAULT_BLACKBODY = 4;
  public static final int DEFAULT_NEON = 5;
  public static final int DEFAULT_WINTER = 6;
  public static final int DEFAULT_SUMMER = 7;


  private class GradientNode {
    float  location; // placement of this node within the gradient - 0.0 to 1.0
    color  col; // colour of the node
    GradientNode(float l, color c) {
      location = l;
      col = c;
    }
  }

  private GradientNode[] nodes;

  Gradient() {
    // create a black to white gradient on construction
    makeDefaultGradient(DEFAULT_BLACK_TO_WHITE);
  }

  void addNode(float location, color col) {
    // add a node to the gradient, specifying the location (0.0-1.0) and the colour
    // NOTE HACK colours should be added in ascending order of location - needs to resort the nodes on adding

    if ( nodes.length == 0 ) {
      nodes = new GradientNode[1];
      nodes[0] = new GradientNode(location, col);
    }
    else {
      GradientNode[] newNodes = new GradientNode[nodes.length+1];
      System.arraycopy(nodes,0, newNodes,0, nodes.length);
      newNodes[nodes.length] = new GradientNode(location, col);
      nodes = newNodes;
    }
  }

  void clear() {
    nodes = new GradientNode[0];
  }

  void makeDefaultGradient(int defaultGradient) {

    clear();

    switch ( defaultGradient ) {
      case DEFAULT_BLACK_TO_WHITE:
        addNode ( 0, 0x000000 );
        addNode ( 1, 0xFFFFFF );  
        break;
      case DEFAULT_RANDOM:
        makeRandomGradient(4);
        break;
      case DEFAULT_SPECTRUM:
        addNode ( 0,    0xFF0000 );
        addNode ( 0.25, 0xFFFF00 );
        addNode ( 0.5,  0x00FF00 );
        addNode ( 0.75, 0x00FFFF );
        addNode ( 1,    0x0000FF );
        break;
      case DEFAULT_INFRARED:
        addNode ( 0,    0x000000 );
        addNode ( 1.0/6,0x000080 );
        addNode ( 2.0/6,0x800080 );
        addNode ( 3.0/6,0xFF0000 );
        addNode ( 4.0/6,0xFF8000 );
        addNode ( 5.0/6,0xFFFF00 );
        addNode ( 1,    0xFFFFFF );
        break;
      case DEFAULT_BLACKBODY:
        addNode ( 0,    0x000000 );
        addNode ( 1.0/5,0x0040FF );
        addNode ( 2.0/5,0x00C0FF );
        addNode ( 3.0/5,0xFF4000 );
        addNode ( 4.0/5,0xFFC000 );
        addNode ( 1,    0xFFFFFF );
        break;
      case DEFAULT_NEON:
        addNode ( 0,    0x000000 );
        addNode ( 0.25, 0x3333FF );
        addNode ( 0.5,  0x0099FF );
        addNode ( 0.75, 0xE60080 );
        addNode ( 1,    0xFF00FF );
        break;
      case DEFAULT_WINTER:
        addNode ( 0,    0x4C80FF );
        addNode ( 0.5,  0xE6E6E6 );
        addNode ( 1,    0x999999 );
        break;
      case DEFAULT_SUMMER:
        addNode ( 0,    0x334CFF );
        addNode ( 0.25, 0xFF0080 );
        addNode ( 0.5,  0xFF8033 );
        addNode ( 0.75, 0xCC4C00 );
        addNode ( 1,    0xFFCC00 );
        break;      
    }

  }

  void makeRandomGradient(int numColours) {
    float location, locationMin, locationMax;
    int r,g,b;

    clear();

    for ( int n=0; n < numColours; n++ ) {
      if ( n == 0 )
        location = 0.0;
      else if ( n == numColours-1 )
        location = 1.0;
      else {
        locationMin = float(n)/numColours;
        locationMax = float(n+1)/numColours;
        location = random ( locationMin, locationMax );
      }

      r = int(random(2.5)) * 128;
      g = int(random(2.5)) * 128;
      b = int(random(2.5)) * 128;

      addNode ( location, color(r,g,b) );
    }
  }

  color getColour(float location) {
    float bandLocation, bandScale, bandDelta;
    float r,g,b;

    for ( int c=0; c < nodes.length-1; c++ ) {
      if ( location >= nodes[c].location && location <= nodes[c+1].location ) {
         bandScale = nodes[c+1].location - nodes[c].location;
         bandLocation = location - nodes[c].location;
         bandDelta = bandLocation / bandScale;

         r = bandDelta * ( red(nodes[c+1].col) - red(nodes[c].col) ) + red(nodes[c].col);
         g = bandDelta * ( green(nodes[c+1].col) - green(nodes[c].col) ) + green(nodes[c].col);
         b = bandDelta * ( blue(nodes[c+1].col) - blue(nodes[c].col) ) + blue(nodes[c].col);
         return color(r,g,b);
      }
    }
    return color(0,0,0);
  }

  color[] makeArrayOfColours(int numColours) {

    float location, bandLocation, bandScale, bandDelta;
    color[] cols = new color[numColours];

    for ( int i=0; i < numColours; i++ ) {
      location = float(i) / (numColours-1);
      cols[i] = getColour(location);
    }

    return cols;
  }

}

tlecoz

Hello !

if I was you, I would divide my 3D object into some 2D-slices (32 slices composed by 32x32 areas for example)

I'll use this wonderfull algorithm http://wonderfl.net/c/7EAn

and I'll apply it on each of my slice. Then "I'll change my point of view" and I'll apply again the algo on each slice (using each slide relative to the previous and the next one)

And I suppose it would work.

But it's a very intensive computation, that's why I used 32x32x32 cube.

Good luck !

prosession

hey, thanks for the answer. The algorithm you have posted is in pure Java code, have you implemented it in processing yet? Can you post the code?

Cheers

tlecoz

Hello ! "The algorithm you have posted is in pure Java code"

Actually, it's not java, it's ActionScript3

I just tryed to translate it in Processing right now, but I get a "null exception" somewhere in the "draw" function and I'm sorry, I'm too busy and have no time to find the bug today.

Here is my actual code, if someone want to try to debug it If nobody do it, I will finish it before the end of the week but don't know when exactly

/**
 * Copyright PaulOllivier ( http//wonderfl.net/user/PaulOllivier )
 * MIT License ( http//www.opensource.org/licenses/mit-license.php )
 * Downloaded from http//wonderfl.net/c/7EAn
 */

int P_NUM;
Particle ppp;
Node[][] grid;
Node nnn;

PImage bmd;

int[] buf;
float mx0;
float my0;


void setup(){
  size(512,512,P2D); 
  bmd = new PImage(512,512);
  bmd.loadPixels();

  buf = new int[512*512];
  for(int i=0;i<512*512;i++)buf[i] = 0;

  grid = new Node[128][128];


  Node node = null;
  int gx;
  int gy;
  for (gx = 0; gx < 128; gx++)
  {
      grid[gx] = new Node[128];

      for (gy = 0; gy < 128; gy++)
      {
          if (nnn != null) grid[gx][gy] = node = node.next = new Node();
          else grid[gx][gy] = node = nnn = new Node();
      }
  }

  for (gx = 0; gx < (128 - 2); gx++)
  {
      for (gy = 0; gy < (128 - 2); gy++)
      {
          grid[gx][gy].n10 = grid[gx + 1][gy + 0];
          grid[gx][gy].n20 = grid[gx + 2][gy + 0];
          grid[gx][gy].n01 = grid[gx + 0][gy + 1];
          grid[gx][gy].n11 = grid[gx + 1][gy + 1];
          grid[gx][gy].n21 = grid[gx + 2][gy + 1];
          grid[gx][gy].n02 = grid[gx + 0][gy + 2];
          grid[gx][gy].n12 = grid[gx + 1][gy + 2];
          grid[gx][gy].n22 = grid[gx + 2][gy + 2];
      }
  }

  Particle p = ppp = new Particle();

  int i = P_NUM;
  while (i-- != 0)
  {
      gx = parseInt( 0.1 + 0.9 * (i & 127));
      gy = parseInt(0.1 + 0.9 * ((i >> 7) & 127));

      if (gx < 25)
      {
          p.col = color(255,0,0);
          p.gas = false;
          p.mass = 4.0;
      } 
      else if (gx < 50)
      {
          p.col = color(255,255,0);
          p.gas = true;
          p.mass = 3.0;
      } 
      else if (gx < 75)
      {
          p.col = color(0,0,255);
          p.gas = false;
          p.mass = 2.0;
      } 
      else
      {
          p.col = color(127,255,0);
          p.gas = false;
          p.mass = 1.0;
      }

      p.rd = 1.0;
      p.k = 1.0;

      p.x = gx + random(1.0) + 4.0;
      p.y = gy + random(1.0) + 4.0;

      p = p.next = new Particle();
  }

}



void draw(){
  float mx = mouseX * 0.25;
  float my = mouseY * 0.25;
  float mdx = mx - mx0;
  float mdy = my - my0;
  mx0 = mx;
  my0 = my;

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - */

  Node n = nnn;
  while (n != null)
  {
      if (n.active)
      {
          n.m = 0;
          n.d = 0;
          n.gx = 0;
          n.gy = 0;
          n.u = 0;
          n.v = 0;
          n.ax = 0;
          n.ay = 0;
          n.active = false;
      }
      n = n.next;
  }

  Particle p = ppp;
  while (p != null)
  {
      int pcx = parseInt(floor(p.x - 0.5) >> 0);
      int pcy = parseInt(floor(p.y - 0.5) >> 0);

      float d = pcx - p.x;
      p.x0 = 0.5 * d * d + 1.5 * d + 1.125;
      p.gx0 = d + 1.5;

      d += 1.0;
      p.x1 = -d * d + 0.75;
      p.gx1 = -2.0 * d;

      d += 1.0;
      p.x2 = 0.5 * d * d - 1.5 * d + 1.125;
      p.gx2 = d - 1.5;

      d = pcy - p.y;
      p.y0 = 0.5 * d * d + 1.5 * d + 1.125;
      p.gy0 = d + 1.5;

      d += 1.0;
      p.y1 = -d * d + 0.75;
      p.gy1 = -2.0 * d;

      d += 1.0;
      p.y2 = 0.5 * d * d - 1.5 * d + 1.125;
      p.gy2 = d - 1.5;

      /* - - - - - - - - - - - - - - - - - - - - - - - - - - - */

      n = p.node = grid[pcx][pcy];

      Node n00 = n;
      Node n01 = n.n01;
      Node n02 = n.n02;
      Node n10 = n.n10;
      Node n11 = n.n11;
      Node n12 = n.n12;
      Node n20 = n.n20;
      Node n21 = n.n21;
      Node n22 = n.n22;

      n00.active = true;
      n01.active = true;
      n02.active = true;
      n10.active = true;
      n11.active = true;
      n12.active = true;
      n20.active = true;
      n21.active = true;
      n22.active = true;

      p.p00 = p.x0 * p.y0;
      p.p10 = p.x1 * p.y0;
      p.p20 = p.x2 * p.y0;
      p.p01 = p.x0 * p.y1;
      p.p11 = p.x1 * p.y1;
      p.p21 = p.x2 * p.y1;
      p.p02 = p.x0 * p.y2;
      p.p12 = p.x1 * p.y2;
      p.p22 = p.x2 * p.y2;

      float pm = p.mass;

      n00.m += p.p00 * pm;
      n00.d += p.p00;
      n00.gx += p.gx0 * p.y0;
      n00.gy += p.x0 * p.gy0;

      n10.m += p.p10 * pm;
      n10.d += p.p10;
      n10.gx += p.gx1 * p.y0;
      n10.gy += p.x1 * p.gy0;

      n20.m += p.p20 * pm;
      n20.d += p.p20;
      n20.gx += p.gx2 * p.y0;
      n20.gy += p.x2 * p.gy0;

      n01.m += p.p01 * pm;
      n01.d += p.p01;
      n01.gx += p.gx0 * p.y1;
      n01.gy += p.x0 * p.gy1;

      n11.m += p.p11 * pm;
      n11.d += p.p11;
      n11.gx += p.gx1 * p.y1;
      n11.gy += p.x1 * p.gy1;

      n21.m += p.p21 * pm;
      n21.d += p.p21;
      n21.gx += p.gx2 * p.y1;
      n21.gy += p.x2 * p.gy1;

      n02.m += p.p02 * pm;
      n02.d += p.p02;
      n02.gx += p.gx0 * p.y2;
      n02.gy += p.x0 * p.gy2;

      n12.m += p.p12 * pm;
      n12.d += p.p12;
      n12.gx += p.gx1 * p.y2;
      n12.gy += p.x1 * p.gy2;

      n22.m += p.p22 * pm;
      n22.d += p.p22;
      n22.gx += p.gx2 * p.y2;
      n22.gy += p.x2 * p.gy2;

      /* - - - - - - - - - - - - - - - - - - - - - - - - - - - */

      float pressure;

      if (p.gas == false)
      {
          d = n10.d - n00.d;

          float C20 = 3.0 * d - n10.gx - 2.0 * n00.gx;
          float C30 = -2.0 * d + n10.gx + n00.gx;

          d = n01.d - n00.d;

          float C02 = 3.0 * d - n01.gy - 2.0 * n00.gy;
          float C03 = -2.0 * d + n01.gy + n00.gy;

          float csum1 = n00.d + n00.gy + C02 + C03;
          float csum2 = n00.d + n00.gx + C20 + C30;

          float C21 = 3.0 * n11.d - 2.0 * n01.gx - n11.gx - 3.0 * csum1 - C20;
          float C31 = -2.0 * n11.d + n01.gx + n11.gx + 2.0 * csum1 - C30;
          float C12 = 3.0 * n11.d - 2.0 * n10.gy - n11.gy - 3.0 * csum2 - C02;
          float C13 = -2.0 * n11.d + n10.gy + n11.gy + 2.0 * csum2 - C03;
          float C11 = n01.gx - C13 - C12 - n00.gx;

          float u = p.x - pcx;
          float u2 = u * u;
          float u3 = u * u2;

          float v = p.y - pcy;
          float v2 = v * v;
          float v3 = v * v2;

          float density = n00.d + n00.gx * u + n00.gy * v + C20 * u2 + C02 * v2 + C30 * u3 + C03 * v3 + C21 * u2 * v + C31 * u3 * v + C12 * u * v2 + C13 * u * v3 + C11 * u * v;

          pressure = p.k * (density - p.rd);
          if (pressure > 2.0) pressure = 2.0;

          pressure *= pm;
      }
      else
      {
          pressure = 1.5;
      }

      /* - - - - - - - - - - - - - - - - - - - - - - - - - - - */

      float fx = 0.0;
      float fy = 0.0;

      if (p.x < 8) fx = pm * (8 - p.x);
      else if (p.x > 128 - 8) fx = pm * (128 - 8 - p.x);

      if (p.y < 8) fy = pm * (8 - p.y);
      else if (p.y > 128 - 8) fy = pm * (128 - 8 - p.y);

      /* - - - - - - - - - - - - - - - - - - - - - - - - - - - */

      float vx = (p.x - mx);
      float vy = (p.y - my);
      if (vx < 0) vx = -vx;
      if (vy < 0) vy = -vy;

      if ((vx < 10.0) && (vy < 10.0))
      {
          float weight = p.mass * (1.0 - (vx / 10.0)) * (1.0 - (vy / 10.0));
          fx += weight * (mdx - p.u);
          fy += weight * (mdy - p.v);
      }

      /* - - - - - - - - - - - - - - - - - - - - - - - - - - - */

      float phi = p.x0 * p.y0;
      n00.ax += -((p.gx0 * p.y0) * pressure) + fx * phi;
      n00.ay += -((p.x0 * p.gy0) * pressure) + fy * phi;

      phi = (p.x0 * p.y1);
      n01.ax += -((p.gx0 * p.y1) * pressure) + fx * phi;
      n01.ay += -((p.x0 * p.gy1) * pressure) + fy * phi;

      phi = (p.x0 * p.y2);
      n02.ax += -((p.gx0 * p.y2) * pressure) + fx * phi;
      n02.ay += -((p.x0 * p.gy2) * pressure) + fy * phi;

      phi = (p.x1 * p.y0);
      n10.ax += -((p.gx1 * p.y0) * pressure) + fx * phi;
      n10.ay += -((p.x1 * p.gy0) * pressure) + fy * phi;

      phi = (p.x1 * p.y1);
      n11.ax += -((p.gx1 * p.y1) * pressure) + fx * phi;
      n11.ay += -((p.x1 * p.gy1) * pressure) + fy * phi;

      phi = (p.x1 * p.y2);
      n12.ax += -((p.gx1 * p.y2) * pressure) + fx * phi;
      n12.ay += -((p.x1 * p.gy2) * pressure) + fy * phi;

      phi = (p.x2 * p.y0);
      n20.ax += -((p.gx2 * p.y0) * pressure) + fx * phi;
      n20.ay += -((p.x2 * p.gy0) * pressure) + fy * phi;

      phi = (p.x2 * p.y1);
      n21.ax += -((p.gx2 * p.y1) * pressure) + fx * phi;
      n21.ay += -((p.x2 * p.gy1) * pressure) + fy * phi;

      phi = (p.x2 * p.y2);
      n22.ax += -((p.gx2 * p.y2) * pressure) + fx * phi;
      n22.ay += -((p.x2 * p.gy2) * pressure) + fy * phi;

      p = p.next;
  }

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

  n = nnn;
  while (n != null)
  {
      if (n.active == true)
      {
          n.ax /= n.m;
          n.ay /= n.m;
          n.ay += 0.03;
      }
      n = n.next;
  }

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

  p = ppp;
  while (p != null)
  {
      n = p.node;

      Node n00 = n;
      Node n01 = n.n01;
      Node n02 = n.n02;
      Node n10 = n.n10;
      Node n11 = n.n11;
      Node n12 = n.n12;
      Node n20 = n.n20;
      Node n21 = n.n21;
      Node n22 = n.n22;

      p.u += (p.p00 * n00.ax + p.p10 * n10.ax + p.p20 * n20.ax + p.p01 * n01.ax + p.p11 * n11.ax + p.p21 * n21.ax + p.p02 * n02.ax + p.p12 * n12.ax + p.p22 * n22.ax);
      p.v += (p.p00 * n00.ay + p.p10 * n10.ay + p.p20 * n20.ay + p.p01 * n01.ay + p.p11 * n11.ay + p.p21 * n21.ay + p.p02 * n02.ay + p.p12 * n12.ay + p.p22 * n22.ay);

      float mu = p.mass * p.u;
      float mv = p.mass * p.v;

      n00.u += p.p00 * mu;
      n00.v += p.p00 * mv;
      n10.u += p.p10 * mu;
      n10.v += p.p10 * mv;
      n20.u += p.p20 * mu;
      n20.v += p.p20 * mv;
      n01.u += p.p01 * mu;
      n01.v += p.p01 * mv;
      n11.u += p.p11 * mu;
      n11.v += p.p11 * mv;
      n21.u += p.p21 * mu;
      n21.v += p.p21 * mv;
      n02.u += p.p02 * mu;
      n02.v += p.p02 * mv;
      n12.u += p.p12 * mu;
      n12.v += p.p12 * mv;
      n22.u += p.p22 * mu;
      n22.v += p.p22 * mv;

      p = p.next;
  }

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

  n = nnn;
  while (n != null)
  {
      if (n.active == true)
      {
          n.u /= n.m;
          n.v /= n.m;
      }
      n = n.next;
  }

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - */


  int pos = 512*512;
  while (pos-->-1)  bmd.pixels[pos] = color(32,18,14);//0x1F1814;

  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

  p = ppp;
  while (p != null)
  {
      n = p.node;

      float gu = p.p00 * n.u + p.p10 * n.n10.u + p.p20 * n.n20.u + p.p01 * n.n01.u + p.p11 * n.n11.u + p.p21 * n.n21.u + p.p02 * n.n02.u + p.p12 * n.n12.u + p.p22 * n.n22.u;
      float gv = p.p00 * n.v + p.p10 * n.n10.v + p.p20 * n.n20.v + p.p01 * n.n01.v + p.p11 * n.n11.v + p.p21 * n.n21.v + p.p02 * n.n02.v + p.p12 * n.n12.v + p.p22 * n.n22.v;

      p.x += gu;
      p.y += gv;
      p.u += gu - p.u;
      p.v += gv - p.v;

      if (p.x < 2)
      {
          p.x = 2.2;
          p.u = 0.0;
      }
      else if (p.x > (128 - 2))
      {
          p.x = 128 - 2.2;
          p.u = 0.0;
      }

      if (p.y < 2)
      {
          p.y = 2.2;
          p.v = 0.0;
      }
      else if (p.y > (128 - 2))
      {
          p.y = 128 - 2.2;
          p.v = 0.0;
      }

      /* - - - - - - - - - - - - - - - - - - - - - - - - - - - */

      int c = p.col;

      int px0 = parseInt(p.x * 4) >> 0;
      int py0 = parseInt(p.y * 4) >> 0;
      int px1 = parseInt((p.x - p.u) * 4) >> 0;
      int py1 = parseInt((p.y - p.v) * 4) >> 0;

      bmd.pixels[px0 + (py0 << 9)] = c;
      bmd.pixels[((px0 + px1) >> 1) + (((py0 + py1) >> 1) << 9)] = c;
      bmd.pixels[px1 + (py1 << 9)] = c;

      p = p.next;
  }

  image(bmd,0,0,512,512);


}



class Node
{
    Node next;
    Boolean active;
    float m = 0;
    float d = 0;
    float gx = 0;
    float gy = 0;
    float u = 0;
    float v = 0;
    float ax = 0;
    float ay = 0;
    Node n10;
    Node n20;
    Node n01;
    Node n11;
    Node n21;
    Node n02;
    Node n12;
    Node n22;

    Node()
    {
    }
}

class Particle
{
      Particle next  ;
      Node node  ;
      int col   = color(127,127,127);//0x808080
      float mass   = 1.0;
      float rd   = 1.0;
      float k   = 1.0;
      Boolean gas;
      float x   = 0.0;
      float y   = 0.0;
      float u   = 0.0;
      float v   = 0.0;
      float gx0  ;
      float gx1  ;
      float gx2  ;
      float gy0  ;
      float gy1  ;
      float gy2  ;
      float x0  ;
      float x1  ;
      float x2  ;
      float y0  ;
      float y1  ;
      float y2  ;
      float p00  ;
      float p10  ;
      float p20  ;
      float p01  ;
      float p11  ;
      float p21  ;
      float p02  ;
      float p12  ;
      float p22  ;

    Particle()
    {
    }
}

prosession

hi, thanks for your work so far! I'll try to figure it out by myself but I'm also excited about your result.

FrankOcean

guy what type of Algo structure are these called ?

    if ( b==1 )  x [0  ][i  ] = -x [1  ][i  ];  else  x [0  ][i  ] = x [1  ][i  ];
    if ( b==1 )  x [n+1][i  ] = -x [n  ][i  ];  else  x [n+1][i  ] = x [n  ][i  ];
    if ( b==2 )  x [i  ][0  ] = -x [i  ][1  ];  else  x [i  ][0  ] = x [i  ][1  ];
    if ( b==2 )  x [i  ][n+1] = -x [i  ][n  ];  else  x [i  ][n+1] = x [i  ][n  ];

or this one?

if ( b==1 )  x [0  ][i  ] = -x [1  ][i  ];  else  x [0  ][i  ] = x [1  ][i  ];
if ( b==1 )  x [n+1][i  ] = -x [n  ][i  ];  else  x [n+1][i  ] = x [n  ][i  ];
if ( b==2 )  x [i  ][0  ] = -x [i  ][1  ];  else  x [i  ][0  ] = x [i  ][1  ];
if ( b==2 )  x [i  ][n+1] = -x [i  ][n  ];  else  x [i  ][n+1] = x [i  ][n  ];

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

i never seen code this neat, AND DAMM thats a long ass code!