176 lines
5.4 KiB
C
176 lines
5.4 KiB
C
void advect64x64(fix* dye, fix* vx, fix* vy, fix* outx, fix* outy, fix* outz) {
|
|
for (int j = 1; j < H-1; j++) {
|
|
for (int i = 1; i < W-1; i++) {
|
|
int id = ID(i, j);
|
|
|
|
fix x = FIX(i) - fmul(dtRdx, vx[id]);
|
|
fix y = FIX(j) - fmul(dtRdx, vy[id]);
|
|
|
|
if (x < 0) x = 0;
|
|
else if (x > MAX_X) x = MAX_X;
|
|
if (y < 0) y = 0;
|
|
else if (y > MAX_Y) y = MAX_Y;
|
|
|
|
int x1 = UNFIX(x);
|
|
int y1 = UNFIX(y);
|
|
int x2 = x1 + 1;
|
|
int y2 = y1 + 1;
|
|
|
|
int tl = ID(x1, y1);
|
|
int tr = ID(x2, y1);
|
|
int bl = ID(x1, y2);
|
|
int br = ID(x2, y2);
|
|
|
|
fix s1 = ffrac(x);
|
|
fix s0 = ONE - s1;
|
|
fix t1 = ffrac(y);
|
|
fix t0 = ONE - t1;
|
|
|
|
outx[id] = fmul(s0, fmul(t0, vx[tl]) + fmul(t1, vx[bl])) +
|
|
fmul(s1, fmul(t0, vx[tr]) + fmul(t1, vx[br]));
|
|
|
|
outy[id] = fmul(s0, fmul(t0, vy[tl]) + fmul(t1, vy[bl])) +
|
|
fmul(s1, fmul(t0, vy[tr]) + fmul(t1, vy[br]));
|
|
|
|
outz[id] = fmul(s0, fmul(t0, dye[tl]) + fmul(t1, dye[bl])) +
|
|
fmul(s1, fmul(t0, dye[tr]) + fmul(t1, dye[br]));
|
|
}
|
|
}
|
|
}
|
|
|
|
void fluid_step64x64() {
|
|
// Advect the velocity field and the dye (27.7ms)
|
|
advect64x64(dye, vx, vy, tx, ty, tz);
|
|
// Swap the pointers
|
|
swap(vx, tx);
|
|
swap(vy, ty);
|
|
swap(dye, tz);
|
|
|
|
// Velocity boundary conditions (0.1ms)
|
|
for (int i = 0; i < W; i++) {
|
|
vy[ID(i, 0)] = -vy[ID(i, 1)];
|
|
vy[ID(i, H-1)] = -vy[ID(i, H-2)];
|
|
vx[ID(i, 0)] = vx[ID(i, 1)];
|
|
vx[ID(i, H-1)] = vx[ID(i, H-2)];
|
|
}
|
|
for (int j = 0; j < H; j++) {
|
|
vx[ID(0, j)] = -vx[ID(1, j)];
|
|
vx[ID(W-1, j)] = -vx[ID(W-2, j)];
|
|
vy[ID(0, j)] = vy[ID(1, j)];
|
|
vy[ID(W-1, j)] = vy[ID(W-2, j)];
|
|
}
|
|
|
|
// Calculate divergence & first pressure iteration (3.6ms)
|
|
for (int j = 1; j < H-1; j++) {
|
|
for (int i = 1; i < W-1; i++) {
|
|
div[ID(i, j)] = fmul(alphaHalfRdx, vx[ID(i+1,j)] - vx[ID(i-1, j)] + vy[ID(i, j+1)] - vy[ID(i, j-1)]);
|
|
p[ID(i, j)] = div[ID(i, j)] >> 2;
|
|
}
|
|
}
|
|
|
|
// Initial Pressure boundary conditions (0ms)
|
|
for (int i = 0; i < W; i++) {
|
|
p[ID(i, 0)] = p[ID(i, 1)];
|
|
p[ID(i, H-1)] = p[ID(i, H-2)];
|
|
}
|
|
for (int j = 0; j < H; j++) {
|
|
p[ID(0, j)] = p[ID(1, j)];
|
|
p[ID(W-1, j)] = p[ID(W-2, j)];
|
|
}
|
|
|
|
// Poisson pressure solver (20.8ms)
|
|
for (int k = 0; k < pressureIterations-1; k++) {
|
|
// Jacobi iteration (20.4ms)
|
|
for (int j = 1; j < H-1; j++) {
|
|
for (int i = 1; i < W-1; i++) {
|
|
int id = ID(i, j);
|
|
p[id] = (div[id] + p[ID(i-1, j)] + p[ID(i+1, j)] + p[ID(i, j-1)] + p[ID(i, j+1)]) >> 2;
|
|
}
|
|
}
|
|
|
|
// Pressure boundary conditions (0.4ms)
|
|
for (int i = 0; i < W; i++) {
|
|
p[ID(i, 0)] = p[ID(i, 1)];
|
|
p[ID(i, H-1)] = p[ID(i, H-2)];
|
|
}
|
|
for (int j = 0; j < H; j++) {
|
|
p[ID(0, j)] = p[ID(1, j)];
|
|
p[ID(W-1, j)] = p[ID(W-2, j)];
|
|
}
|
|
}
|
|
|
|
// Apply Forces
|
|
if (currentView == 0) {
|
|
int lastX = addX;
|
|
int lastY = addY;
|
|
|
|
if (keydown(KEY_LEFT)) { addX -= 1; if (addX < 1) addX = W-2; }
|
|
else if (keydown(KEY_RIGHT)) { addX += 1; if (addX > W-2) addX = 1; }
|
|
if (keydown(KEY_UP)) { addY -= 1; if (addY < 1) addY = H-2; }
|
|
else if (keydown(KEY_DOWN)) { addY += 1; if (addY > H-2) addY = 1; }
|
|
|
|
fix adx = FIX(addX - lastX);
|
|
fix ady = FIX(addY - lastY);
|
|
|
|
if (adx || ady) {
|
|
fix ddx = abs(adx);
|
|
fix ddy = abs(ady);
|
|
fix vdx = adx / 10;
|
|
fix vdy = ady / 10;
|
|
fix ddxy = ddx + ddy;
|
|
|
|
for (int j = addY-3; j < addY+3; j++) {
|
|
for (int i = addX-3; i < addX+3; i++) {
|
|
if (i<=0||j<=0||i>=W-1||j>=H-1) continue;
|
|
int id = ID(i, j);
|
|
|
|
fix dist = radius - FIX((addX-i)*(addX-i)+(addY-j)*(addY-j));
|
|
if (dist < 0) continue;
|
|
fix dyeD = fmul(dist, dyeIntensity);
|
|
fix velD = fmul(dist, velIntensity);
|
|
|
|
dye[id] = dye[id] + fmul(ddxy, dyeD);
|
|
if (dye[id] > maxDyeIntensity) dye[id] = maxDyeIntensity;
|
|
vx[id] = vx[id] + fmul(vdx, velD);
|
|
vy[id] = vy[id] + fmul(vdy, velD);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// Subtract pressure from velocity & Diffuse dye and velocity (3.4ms)
|
|
for (int j = 1; j < H-1; j++) {
|
|
for (int i = 1; i < W-1; i++) {
|
|
int id = ID(i, j);
|
|
vx[id] = fmul(vx[id] - fmul(halfRdx, p[ID(i+1, j)] - p[ID(i-1, j)]), velDiffusion);
|
|
vy[id] = fmul(vy[id] - fmul(halfRdx, p[ID(i, j+1)] - p[ID(i, j-1)]), velDiffusion);
|
|
}
|
|
}
|
|
|
|
// Draw dye to vram (3.1ms)
|
|
uint32_t *vram = gint_vram+4;
|
|
for (int j = 1; j < H-1; j++) {
|
|
uint32_t data = 0;
|
|
|
|
data = (data << 1);
|
|
for(int k = 1; k < 32; k++) {
|
|
int x = 32*0+k;
|
|
int id = ID(x, j);
|
|
dye[id] = fmul(dye[id], dyeDiffusion);
|
|
data = (data << 1) | (dye[id] > saturateTreshold || (dye[id] > emptyTreshold && (x+j)%2));
|
|
}
|
|
vram[0] = data;
|
|
data = 0;
|
|
|
|
for(int k = 0; k < 31; k++) {
|
|
int x = 32*1+k;
|
|
int id = ID(x, j);
|
|
dye[id] = fmul(dye[id], dyeDiffusion);
|
|
data = (data << 1) | (dye[id] > saturateTreshold || (dye[id] > emptyTreshold && (x+j)%2));
|
|
}
|
|
data = (data << 1);
|
|
vram[1] = data;
|
|
|
|
vram+=4;
|
|
}
|
|
} |