- changed splat to adapt to canvas

- compare stencils script
- blit divergence vertex shader

compare_stencils.py

@@ -0,0 +1,106 @@
+
+from sympy import sqrt, symbols, factorial, IndexedBase, Symbol
+
+p = IndexedBase("p")
+h = Symbol("h")
+
+def taylor(a, b, c):
+    n = 4
+    derivatives = [0] * n
+    steps = [0] * n
+
+    result = p
+
+    for order in range(1, n+1):
+        for i in range(0, 3**order):
+            for j in range(0, order):
+                derivatives[j] = int((i / (3**(order-1-j))) % 3 + 1)
+                step = [a, b, c][derivatives[j] - 1]
+                steps[j] = step
+
+            total_step = 1
+            symbol_name = ""
+            for j in range(0, order):
+                total_step *= (steps[j] * h)
+                component = derivatives[j] - 1
+                symbol_name += ["x", "y", "z"][component]
+
+            result += (1 / factorial(order)) * total_step * Symbol(symbol_name)
+    return result;
+
+def expand_laplacian(face, corner):
+    # Note: Any linear combination of the face and corner coefficients will
+    # produce a valid discretized Laplacian, except for when normalization == 0
+    center = face * 4 + corner * 4
+    normalization = face + 2 * corner
+
+    expr = (
+        - center * taylor(0, 0, 0)
+
+        + face * taylor( 1,  0,  0)
+        + face * taylor(-1,  0,  0)
+        + face * taylor( 0, -1,  0)
+        + face * taylor( 0,  1,  0)
+
+        +  corner * taylor( 1,  1,  0)
+        +  corner * taylor( 1, -1,  0)
+        +  corner * taylor(-1,  1,  0)
+        +  corner * taylor(-1, -1,  0)
+
+    ) / (normalization*h*h)
+
+    print("Laplacian:")
+    print("\tCenter:        {:+d}".format(-center))
+    print("\tFace:          {:+d}".format(face))
+    print("\tCorner:        {:+d}".format(corner))
+    print("\tNormalization: {:+d}".format(normalization))
+    print("\tExpansion:     {}".format(expr.simplify()))
+    print()
+
+    return expr
+
+#a = expand_laplacian(1, 0)
+#b = expand_laplacian(0, 1)
+#c = expand_laplacian(1, 1)
+#d = expand_laplacian(2, 1)
+#e = expand_laplacian(1, 2)
+
+expr_colocated = (
+    -4*taylor(0, 0, 0)
+    +1*taylor(-2, 0, 0)
+    +1*taylor(+2, 0, 0)
+    +1*taylor(0, -2, 0)
+    +1*taylor(0, +2, 0)
+) / (4*h*h)
+
+expr_wide = (
+    -8*taylor(0, 0, 0)
+    +1*taylor(-1, 0, 0)
+    +1*taylor(+1, 0, 0)
+    +1*taylor(0, -1, 0)
+    +1*taylor(0, +1, 0)
+    +1*taylor(-1,-1, 0)
+    +1*taylor(+1,+1, 0)
+    +1*taylor(+1, -1, 0)
+    +1*taylor(-1, +1, 0)
+) / (3*h*h)
+
+expr_vertex = (
+    -4*taylor(0, 0, 0)
+    +1*taylor(-1,-1, 0)
+    +1*taylor(+1,+1, 0)
+    +1*taylor(+1, -1, 0)
+    +1*taylor(-1, +1, 0)
+) / (2*h*h)
+
+expr_mac = (
+    -4*taylor(0, 0, 0)
+    +1*taylor(-1, 0, 0)
+    +1*taylor(+1, 0, 0)
+    +1*taylor(0, -1, 0)
+    +1*taylor(0, +1, 0)
+) / (h*h)
+
+print((expr_colocated - expr_mac).simplify())
+print((expr_colocated - expr_wide).simplify())
+print((expr_colocated - expr_vertex).simplify())

doc/notes.txt

@@ -0,0 +1,56 @@
+
+Colocated multigrid
+
+    +---+---+---+---+---+---+---+---+---+
+ 8  | x | x | x | x | x | x | x | x | x |
+    +---+---+---+---+---+---+---+---+---+
+ 7  | 0 |   |   |   |   |   |   |   | x |
+    +---+---+---+---+---+---+---+---+---+
+ 6  | 0 |   |   |   |   |   |   |   | x |
+    +---+---+---+---+---+---+---+---+---+
+ 5  | 0 |   |   |   |   |   |   |   | x |
+    +---+---+---+---+---+---+---+---+---+
+ 4  | 0 |   |   |   |   |   |   |   | x |
+    +---+---+---+---+---+---+---+---+---+
+ 3  | 0 |   |   |   |   |   |   |   | x |
+    +---+---+---+---+---+---+---+---+---+
+ 2  | 0 |   |   |   |   |   |   |   | x |
+    +---+---+---+---+---+---+---+---+---+
+ 1  | 0 |   |   |   |   |   |   |   | x |
+    +---+---+---+---+---+---+---+---+---+
+ 0  | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | x |
+    +---+---+---+---+---+---+---+---+---+
+      0   1   2   3   4   5   6   7   8
+
+- Pressure is defined on all cells:
+	- The 0 cells are stored but not used (we return p = 0 on these cells)
+	- The x cells are not stored (we implicitly return p = 0 out of bounds)
+
+- Velocity is only defined in the inner cells, we always return 0 on the borders
+
+The Coarse grid is as follows:
+
+  +-------+-------+-------+-------+-------+
+  |       |       |       |       |       |
+  |   x   |   x   |   x   |   x   |   x   |
+  |       |       |       |       |       |
+  +-------+-------+-------+-------+-------+
+  |       |       |       |       |       |
+  |   0   |       |       |       |   x   |
+  |       |       |       |       |       |
+  +-------+-------+-------+-------+-------+
+  |       |       |       |       |       |
+  |   0   |       |       |       |   x   |
+  |       |       |       |       |       |
+  +-------+-------+-------+-------+-------+
+  |       |       |       |       |       |
+  |   0   |       |       |       |   x   |
+  |       |       |       |       |       |
+  +-------+-------+-------+-------+-------+
+  |       |       |       |       |       |
+  |   0   |   0   |   0   |   0   |   x   |
+  |       |       |       |       |       |
+  +-------+-------+-------+-------+-------+
+
+So, the size of the total grid is (2^N+1)*(2^N+1), but we store only (2^N)*(2^N), and explicitly
+set the pressure to 0 on the lower and left borders.

loader.js

@@ -808,6 +808,7 @@ function GLBindImports(env)
 
 	env.glUniform1f = function(loc, v0) { GLctx.uniform1f(GLuniforms[loc], v0); };
 	env.glUniform1i = function(loc, v0) { GLctx.uniform1i(GLuniforms[loc], v0); };
+	env.glUniform2i = function(loc, v0, v1) { GLctx.uniform2i(GLuniforms[loc], v0, v1); };
 	env.glUniform2f = function(loc, v0, v1) { GLctx.uniform2f(GLuniforms[loc], v0, v1); };
 	env.glUniform3f = function(loc, v0, v1, v2) { GLctx.uniform3f(GLuniforms[loc], v0, v1, v2); };
 

notes.txt

@@ -1,11 +0,0 @@
-On macOS, homebrew version of llvm must be used instead of xcode's one:
-
-> brew install llvm
-
-Have to run a local http server to serve files:
-
-> python3 -m http.server
-
-
-brew llvm toolchain seems to be broken and does not find open gl header files, and trying to include them from the framework path in /applications/xcode/ etc creates a bunch of other errors, so the workaround was to download mesa openGL headers, and tweak glext.h so that it includes stdint.h instead of inttypes.h (because brew's llvm can't seem to find this one). This was annoying as fuck.
-Zig cc seem to have everything needed but doesn't recognize --export-all and such, so I could'nt find a way to use it.

src/main.c

@@ -16,6 +16,8 @@
 
 #include"blit_vertex_shader.h"
 #include"blit_fragment_shader.h"
+#include"blit_div_vertex_shader.h"
+#include"blit_div_fragment_shader.h"
 #include"common_vertex_shader.h"
 #include"advect_shader.h"
 #include"divergence_shader.h"
@@ -24,6 +26,7 @@
 #include"multigrid_correct_shader.h"
 #include"subtract_pressure_shader.h"
 #include"splat_shader.h"
+#include"blit_residue_fragment_shader.h"
 
 void* memcpy(void* dst, void* src, unsigned long n)
 {
@@ -50,7 +53,7 @@ float fabs(float x)
 	return(x>0 ? x : -x);
 }
 
-#define PI 3.14159
+#define PI 3.14159265359
 #define EULER 2.71828
 #define SQRT2 1.41421
 
@@ -133,6 +136,16 @@ typedef struct jacobi_program
 
 } jacobi_program;
 
+typedef struct blit_residue_program
+{
+	GLuint prog;
+
+	GLint pos;
+	GLint mvp;
+	GLint xTex;
+	GLint bTex;
+} blit_residue_program;
+
 typedef struct multigrid_restrict_residual_program
 {
 	GLuint prog;
@@ -169,6 +182,7 @@ typedef struct blit_program
 
 	GLint pos;
 	GLint mvp;
+	GLint gridSize;
 	GLint tex;
 } blit_program;
 
@@ -183,6 +197,7 @@ typedef struct splat_program
 	GLint radius;
 	GLint additive;
 	GLint blending;
+	GLint randomize;
 
 } splat_program;
 
@@ -201,7 +216,8 @@ multigrid_correct_program multigridCorrectProgram;
 subtract_program subtractProgram;
 splat_program splatProgram;
 blit_program blitProgram;
-
+blit_program blitDivProgram;
+blit_residue_program blitResidueProgram;
 
 frame_buffer colorBuffer;
 frame_buffer velocityBuffer;
@@ -302,6 +318,7 @@ void init_splat(splat_program* program)
 	program->radius = glGetUniformLocation(program->prog, "radius");
 	program->additive = glGetUniformLocation(program->prog, "additive");
 	program->blending = glGetUniformLocation(program->prog, "blending");
+	program->randomize = glGetUniformLocation(program->prog, "randomize");
 }
 
 void init_blit(blit_program* program)
@@ -311,8 +328,29 @@ void init_blit(blit_program* program)
 	program->pos = glGetAttribLocation(program->prog, "pos");
 	program->mvp = glGetUniformLocation(program->prog, "mvp");
 	program->tex = glGetUniformLocation(program->prog, "tex");
+	program->gridSize = glGetUniformLocation(program->prog, "gridSize");
+}
+
+void init_blit_div(blit_program* program)
+{
+	console_log_fmt("compiling blit div...");
+	program->prog = compile_shader(blit_div_vertex_shader_src, blit_div_fragment_shader_src);
+	program->pos = glGetAttribLocation(program->prog, "pos");
+	program->mvp = glGetUniformLocation(program->prog, "mvp");
+	program->tex = glGetUniformLocation(program->prog, "tex");
 }
 
+void init_blit_residue(blit_residue_program* program)
+{
+	console_log_fmt("compiling blit residue...");
+	program->prog = compile_shader(blit_div_vertex_shader_src, blit_residue_fragment_shader_src);
+	program->pos = glGetAttribLocation(program->prog, "pos");
+	program->mvp = glGetUniformLocation(program->prog, "mvp");
+	program->xTex = glGetUniformLocation(program->prog, "xTex");
+	program->bTex = glGetUniformLocation(program->prog, "bTex");
+}
+
+
 GLuint create_texture(int width, int height, GLenum internalFormat, GLenum format, GLenum type, char* initData)
 {
 	GLuint texture;
@@ -323,7 +361,6 @@ GLuint create_texture(int width, int height, GLenum internalFormat, GLenum forma
 	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR);
 	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
 	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
-
 	return(texture);
 }
 
@@ -334,7 +371,6 @@ GLuint create_fbo(GLuint texture)
 	glBindFramebuffer(GL_FRAMEBUFFER, fbo);
 	glBindTexture(GL_TEXTURE_2D, texture);
 	glFramebufferTexture2D(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT0, GL_TEXTURE_2D, texture, 0);
-	glClear(GL_COLOR_BUFFER_BIT);
 	return(fbo);
 }
 
@@ -373,10 +409,7 @@ void frame_buffer_swap(frame_buffer* buffer)
 //----------------------------------------------------------------
 //NOTE(martin): entry point
 //----------------------------------------------------------------
-/*
-#define texWidth 512
-#define texHeight 512
-*/
+
 #define texWidth (256)
 #define texHeight (256)
 
@@ -397,8 +430,11 @@ void reset_texture(GLuint texture, float width, float height, char* initData)
 	glTexImage2D(GL_TEXTURE_2D, 0, TEX_INTERNAL_FORMAT, width, height, 0, TEX_FORMAT, TEX_TYPE, initData);
 }
 
+//static bool resetCmd = false;
+
 void reset()
 {
+//	resetCmd = true;
 	console_log_fmt("reset");
 
 	reset_texture(colorBuffer.textures[0], texWidth, texHeight, (char*)colorInitData);
@@ -476,6 +512,7 @@ void init_velocity_vortex()
 			float y = 2*i/(float)texWidth - 1;
 			velocityInitData[i][j][0] = sinf(2*PI*y);
 			velocityInitData[i][j][1] = sinf(2*PI*x);
+
 		}
 	}
 }
@@ -490,22 +527,27 @@ int init(float canvasSize)
 
 	WAJS_SetupCanvas(canvasSize, canvasSize);
 
+//	init_color_checker();
+//	init_velocity_vortex();
+
 	// init programs
 	init_advect(&advectProgram);
 	init_div(&divProgram);
 	init_jacobi(&jacobiProgram);
 	init_multigrid_restrict_residual(&multigridRestrictResidualProgram);
 	init_multigrid_correct(&multigridCorrectProgram);
+	init_blit_residue(&blitResidueProgram);
 
 	init_subtract(&subtractProgram);
 	init_splat(&splatProgram);
 	init_blit(&blitProgram);
+	init_blit_div(&blitDivProgram);
 
 	// init frame buffers
 	console_log_fmt("create color buffer");
-	init_frame_buffer(&colorBuffer, texWidth, texHeight, TEX_INTERNAL_FORMAT, TEX_FORMAT, TEX_TYPE, 0);
+	init_frame_buffer(&colorBuffer, texWidth, texHeight, TEX_INTERNAL_FORMAT, TEX_FORMAT, TEX_TYPE, (char*)colorInitData);
 	console_log_fmt("create velocity buffer");
-	init_frame_buffer(&velocityBuffer, texWidth, texHeight, TEX_INTERNAL_FORMAT, TEX_FORMAT, TEX_TYPE, 0);
+	init_frame_buffer(&velocityBuffer, texWidth, texHeight, TEX_INTERNAL_FORMAT, TEX_FORMAT, TEX_TYPE, (char*)velocityInitData);
 
 	int gridFactor = 1;
 	for(int i=0; i<MULTIGRID_COUNT; i++)
@@ -547,10 +589,19 @@ int init(float canvasSize)
 	return(0);
 }
 
-void apply_splat(float splatPosX, float splatPosY, float splatVelX, float splatVelY, float r, float g, float b)
+void apply_splat(float splatPosX, float splatPosY, float radius, float splatVelX, float splatVelY, float r, float g, float b, bool randomize)
 {
 	glUseProgram(splatProgram.prog);
 
+	if(randomize)
+	{
+		glUniform1f(splatProgram.randomize, 1.);
+	}
+	else
+	{
+		glUniform1f(splatProgram.randomize, 0.);
+	}
+
 	// force
 	glBindFramebuffer(GL_FRAMEBUFFER, velocityBuffer.fbos[1]);
 
@@ -563,7 +614,7 @@ void apply_splat(float splatPosX, float splatPosY, float splatVelX, float splatV
 	glUniform1f(splatProgram.additive, 1);
 	glUniform1f(splatProgram.blending, 0);
 
-	glUniform1f(splatProgram.radius, 0.01);
+	glUniform1f(splatProgram.radius, radius);
 
 	glDrawArrays(GL_TRIANGLES, 0, 6);
 
@@ -580,7 +631,7 @@ void apply_splat(float splatPosX, float splatPosY, float splatVelX, float splatV
 	glUniform3f(splatProgram.splatColor, r, g, b);
 	glUniform1f(splatProgram.additive, 0);
 	glUniform1f(splatProgram.blending, 1);
-	glUniform1f(splatProgram.radius, 0.01);
+	glUniform1f(splatProgram.radius, radius);
 
 	glDrawArrays(GL_TRIANGLES, 0, 6);
 
@@ -588,7 +639,7 @@ void apply_splat(float splatPosX, float splatPosY, float splatVelX, float splatV
 
 }
 
-void multigrid_smooth(frame_buffer* x, frame_buffer* b, float invGridSize, int iterationCount)
+void jacobi_solve(frame_buffer* x, frame_buffer* b, float invGridSize, int iterationCount)
 {
 	glUseProgram(jacobiProgram.prog);
 
@@ -656,6 +707,40 @@ void multigrid_clear(frame_buffer* error)
 	glClear(GL_COLOR_BUFFER_BIT);
 }
 
+void input_splat(float t)
+{
+	//NOTE: apply force and dye
+	if(mouseInput.down && (mouseInput.deltaX || mouseInput.deltaY))
+	{
+		// account for margin
+		float margin = 32;
+
+		float canvasWidth = canvas_width();
+		float canvasHeight = canvas_height();
+		float offset = margin/texWidth;
+		float ratio = 1 - 2*margin/texWidth;
+
+		float splatPosX = (mouseInput.x/canvasWidth)*ratio + offset;
+		float splatPosY = (1 - mouseInput.y/canvasHeight)*ratio + offset;
+
+		float splatVelX = (10000.*DELTA*mouseInput.deltaX/canvasWidth)*ratio;
+		float splatVelY = (-10000.*DELTA*mouseInput.deltaY/canvasWidth)*ratio;
+
+		float intensity = 100*sqrtf(square(ratio*mouseInput.deltaX/canvasWidth) + square(ratio*mouseInput.deltaY/canvasHeight));
+
+		float r = intensity * (sinf(2*PI*0.1*t) + 1);
+		float g = 0.5*intensity * (cosf(2*PI*0.1/EULER*t + 654) + 1);
+		float b = intensity * (sinf(2*PI*0.1/SQRT2*t + 937) + 1);
+
+		float radius = 0.005;
+
+		apply_splat(splatPosX, splatPosY, radius, splatVelX, splatVelY, r, g, b, false);
+
+		mouseInput.deltaX = 0;
+		mouseInput.deltaY = 0;
+	}
+}
+
 // This function is called by loader.js every frame
 void WAFNDraw()
 {
@@ -686,39 +771,60 @@ void WAFNDraw()
 
 	frame_buffer_swap(&velocityBuffer);
 
-	//NOTE: apply force and dye
-	if(mouseInput.down && (mouseInput.deltaX || mouseInput.deltaY))
-	{
-		float canvasWidth = canvas_width();
-		float canvasHeight = canvas_height();
-
-		float splatPosX = mouseInput.x/canvasWidth;
-		float splatPosY = 1 - mouseInput.y/canvasHeight;
-
-		float splatVelX = 10000.*DELTA*mouseInput.deltaX/canvasWidth;
-		float splatVelY = -10000.*DELTA*mouseInput.deltaY/canvasWidth;
-
-		float intensity = 100*sqrtf(square(mouseInput.deltaX/canvasWidth) + square(mouseInput.deltaY/canvasHeight));
+	/*
+	//DEBUG
+	static bool splatTrig = false;
+	static bool splat = false;
+	static float splatStart = 0;
+	static int splatDir = 0;
 
-		float r = intensity * (sinf(2*PI*0.1*t) + 1);
-		float g = intensity * (cosf(2*PI*0.1/EULER*t + 654) + 1);
-		float b = intensity * (sinf(2*PI*0.1/SQRT2*t + 937) + 1);
+	static int frameCount = 0;
 
-		apply_splat(splatPosX, splatPosY, splatVelX, splatVelY, r, g, b);
+	if(resetCmd)
+	{
+		frameCount = 0;
+		splat = true;
+		splatStart = frameT;
+	}
 
-		mouseInput.deltaX = 0;
-		mouseInput.deltaY = 0;
+	if(splat)
+	{
+		if(frameT - splatStart >= 0.5)
+		{
+			splat = false;
+			splatDir++;
+			splatDir = splatDir % 3;
+		}
+		float dirX = 0;
+		float dirY = 0;
+		if(splatDir == 0)
+		{
+			dirX = 0;
+			dirY = 0.3;
+		}
+		if(splatDir == 1)
+		{
+			dirX = 0.3;
+			dirY = 0;
+		}
+		if(splatDir == 2)
+		{
+			dirX = 0.2121;
+			dirY = 0.2121;
+		}
+		apply_splat(0.5, 0.5, dirX, dirY, 1.5, 1., 0.1, false);
 	}
+	resetCmd = false;
 
-	//DEBUG
-	//*
-	if((int)(frameT*3/2.) % 6 == 0)
+	if(frameCount>20)
 	{
-		float dirX = 0.2*cosf(2*3.141516*frameT*7.26);
-		float dirY = 0.3+0.1*sinf(frameT);
-		apply_splat(0.5, 0.1, dirX, dirY, 1.5, 1., 0.1);
+		return;
 	}
-	//*/
+	frameCount++;
+*/
+
+	input_splat(t);
+
 
 	//NOTE: compute divergence of advected velocity
 	glUseProgram(divProgram.prog);
@@ -733,48 +839,37 @@ void WAFNDraw()
 	frame_buffer_swap(&divBuffer[0]);
 
 	//NOTE: compute pressure
+	glBindFramebuffer(GL_FRAMEBUFFER, pressureBuffer[0].fbos[1]);
+	glClear(GL_COLOR_BUFFER_BIT);
 
 	#if 0
 		multigrid_clear(&pressureBuffer[0]);
-		multigrid_smooth(&pressureBuffer[0], &divBuffer[0], INV_GRID_SIZE, 500);
+		jacobi_solve(&pressureBuffer[0], &divBuffer[0], INV_GRID_SIZE, texWidth*texHeight);
 	#else
-	//*
 		multigrid_clear(&pressureBuffer[0]);
 
-		for(int i=0; i<4; i++)
+		for(int i=0; i<1; i++)
 		{
-			multigrid_smooth(&pressureBuffer[0], &divBuffer[0], INV_GRID_SIZE, 2);
-			//*
+			jacobi_solve(&pressureBuffer[0], &divBuffer[0], INV_GRID_SIZE, 2);
 			multigrid_coarsen_residual(&divBuffer[1], &pressureBuffer[0], &divBuffer[0], INV_GRID_SIZE);
 
 			multigrid_clear(&pressureBuffer[1]);
-			multigrid_smooth(&pressureBuffer[1], &divBuffer[1], 2*INV_GRID_SIZE, 2);
+			jacobi_solve(&pressureBuffer[1], &divBuffer[1], 2*INV_GRID_SIZE, 2);
 			multigrid_coarsen_residual(&divBuffer[2], &pressureBuffer[1], &divBuffer[1], 2*INV_GRID_SIZE);
 
 			multigrid_clear(&pressureBuffer[2]);
-			multigrid_smooth(&pressureBuffer[2], &divBuffer[2], 4*INV_GRID_SIZE, 60);
+			jacobi_solve(&pressureBuffer[2], &divBuffer[2], 4*INV_GRID_SIZE, 30);
 
 			multigrid_prolongate_and_correct(&pressureBuffer[1], &pressureBuffer[2], 2*INV_GRID_SIZE);
-			multigrid_smooth(&pressureBuffer[1], &divBuffer[1], 2*INV_GRID_SIZE, 8);
+			jacobi_solve(&pressureBuffer[1], &divBuffer[1], 2*INV_GRID_SIZE, 8);
 
 			multigrid_prolongate_and_correct(&pressureBuffer[0], &pressureBuffer[1], INV_GRID_SIZE);
-			//*/
-			multigrid_smooth(&pressureBuffer[0], &divBuffer[0], INV_GRID_SIZE, 4);
+			jacobi_solve(&pressureBuffer[0], &divBuffer[0], INV_GRID_SIZE, 4);
 		}
-		/*/
-		multigrid_clear(&pressureBuffer[0]);
-		multigrid_smooth(&pressureBuffer[0], &pressureBuffer[0], &divBuffer[0], INV_GRID_SIZE, 10);
-
-		multigrid_coarsen_residual(&divBuffer[1], &pressureBuffer[0], &divBuffer[0], INV_GRID_SIZE);
+	#endif
 
-		multigrid_clear(&pressureBuffer[1]);
-		multigrid_smooth(&pressureBuffer[1], &pressureBuffer[1], &divBuffer[1], 2*INV_GRID_SIZE, 400);
 
-		multigrid_prolongate_and_correct(&pressureBuffer[0], &pressureBuffer[0], &pressureBuffer[1], INV_GRID_SIZE);
 
-		multigrid_smooth(&pressureBuffer[0], &pressureBuffer[0], &divBuffer[0], INV_GRID_SIZE, 10);
-	//*/
-	#endif
 
 	//NOTE: subtract pressure gradient to advected velocity
 	glUseProgram(subtractProgram.prog);
@@ -808,23 +903,40 @@ void WAFNDraw()
 
 	glUniform1f(advectProgram.delta, DELTA);
 
-	glUniform1f(advectProgram.dissipation, 0.5);
+	glUniform1f(advectProgram.dissipation, 0.001);
 
 	glDrawArrays(GL_TRIANGLES, 0, 6);
 
 	frame_buffer_swap(&colorBuffer);
 
-//*
 	//NOTE: Blit color texture to screen
 
+	//NOTE: blit residue to screen
+	glViewport(0, 0, canvas_width(), canvas_height());
+
 	float displayMatrix[16] = {
 		1/aspectRatio, 0, 0, 0,
 		0, 1, 0, 0,
 		0, 0, 1, 0,
 		0, 0, 0, 1 };
 
-	glViewport(0, 0, canvas_width(), canvas_height());
+/*
+	glUseProgram(blitResidueProgram.prog);
+	glBindFramebuffer(GL_FRAMEBUFFER, 0);
+
+	glActiveTexture(GL_TEXTURE0);
+	glBindTexture(GL_TEXTURE_2D, pressureBuffer[0].textures[0]);
+	glUniform1i(blitResidueProgram.xTex, 0);
+
+	glActiveTexture(GL_TEXTURE1);
+	glBindTexture(GL_TEXTURE_2D, divBuffer[0].textures[0]);
+	glUniform1i(blitResidueProgram.bTex, 1);
 
+	glUniformMatrix4fv(blitResidueProgram.mvp, 1, GL_FALSE, displayMatrix);
+
+	glDrawArrays(GL_TRIANGLES, 0, 6);
+//*/
+//*
 	glUseProgram(blitProgram.prog);
 	glBindFramebuffer(GL_FRAMEBUFFER, 0);
 
@@ -832,21 +944,37 @@ void WAFNDraw()
 	glBindTexture(GL_TEXTURE_2D, colorBuffer.textures[0]);
 	glUniform1i(blitProgram.tex, 0);
 
+	glUniform2i(blitProgram.gridSize, texWidth, texHeight);
+
 	glUniformMatrix4fv(blitProgram.mvp, 1, GL_FALSE, displayMatrix);
 
 	glDrawArrays(GL_TRIANGLES, 0, 6);
 /*/
-	//NOTE: Blit velocity to screen
-	glViewport(0, 0, window_width(), window_height());
-	glUseProgram(blitProgram.prog);
-	glBindFramebuffer(GL_FRAMEBUFFER, 0);
+
+	//NOTE: recompute divergence of (corrected) velocity
+	glUseProgram(divProgram.prog);
+	glBindFramebuffer(GL_FRAMEBUFFER, divBuffer[0].fbos[1]);
 
 	glActiveTexture(GL_TEXTURE0);
 	glBindTexture(GL_TEXTURE_2D, velocityBuffer.textures[0]);
-	glUniform1i(blitProgram.tex, 0);
+	glUniform1i(divProgram.src, 0);
 
-	glUniformMatrix4fv(blitProgram.mvp, 1, GL_FALSE, displayMatrix);
+	glDrawArrays(GL_TRIANGLES, 0, 6);
+
+	frame_buffer_swap(&divBuffer[0]);
+
+	//NOTE: Blit divergence to screen
+	glViewport(0, 0, canvas_width(), canvas_height());
+	glUseProgram(blitDivProgram.prog);
+	glBindFramebuffer(GL_FRAMEBUFFER, 0);
+
+	glActiveTexture(GL_TEXTURE0);
+	glBindTexture(GL_TEXTURE_2D, divBuffer[0].textures[0]);
+	glUniform1i(blitDivProgram.tex, 0);
+
+	glUniformMatrix4fv(blitDivProgram.mvp, 1, GL_FALSE, displayMatrix);
 
 	glDrawArrays(GL_TRIANGLES, 0, 6);
+
 //*/
 }

src/shaders/advect.glsl

@@ -11,9 +11,46 @@ uniform sampler2D velocity;
 uniform float delta;
 uniform float dissipation;
 
+vec2 u(ivec2 coord)
+{
+	return(texelFetch(velocity, coord, 0).xy);
+}
+
+vec4 q(ivec2 coord)
+{
+	if(  coord.x < 0
+	  || coord.x >= textureSize(src, 0).x
+	  || coord.y < 0
+	  || coord.y >= textureSize(src, 0).y)
+	{
+		return(vec4(0.));
+	}
+	return(texelFetch(src, coord, 0));
+}
+
+vec4 bilerpSrc(vec2 pos)
+{
+	vec2 offset = fract(pos);
+
+	ivec2 bl = ivec2(floor(pos));
+
+	ivec2 br = bl + ivec2(1, 0);
+	ivec2 tl = bl + ivec2(0, 1);
+	ivec2 tr = bl + ivec2(1, 1);
+
+	vec4 lerpTop = (1.-offset.x)*q(tl) + offset.x*q(tr);
+	vec4 lerpBottom = (1.-offset.x)*q(bl) + offset.x*q(br);
+	vec4 result = (1.-offset.y)*lerpBottom + offset.y*lerpTop;
+
+	return(result);
+}
+
 void main()
 {
-	vec2 u = texture(velocity, texCoord).xy;
-	vec2 samplePos = texCoord - delta * u;
-	fragColor = texture(src, samplePos) / (1. + dissipation*delta);
+	float texWidth = float(textureSize(velocity, 0).x);
+
+	ivec2 pixelCoord = ivec2(floor(gl_FragCoord.xy));
+
+	vec2 samplePos = vec2(pixelCoord) - texWidth * delta * u(pixelCoord);
+	fragColor = bilerpSrc(samplePos) / (1. + dissipation*delta);
 }

src/shaders/blit_div_fragment.glsl

@@ -0,0 +1,34 @@
+#version 300 es
+
+precision highp float;
+precision highp sampler2D;
+
+in vec2 texCoord;
+out vec4 fragColor;
+
+uniform sampler2D tex;
+
+vec3 color_map(float v)
+{
+	float logv = log(abs(v))/log(10.0);
+	float f = floor(logv + 7.0);
+	float i = floor(4.0*(logv + 7.0 - f));
+
+	if(f < 0.0) return vec3(0.0);
+	if(f < 1.0) return mix(vec3(1.0, 0.0, 0.0), vec3(1.0), i/4.0);
+	if(f < 2.0) return mix(vec3(0.0, 1.0, 0.0), vec3(1.0), i/4.0);
+	if(f < 3.0) return mix(vec3(0.0, 0.0, 1.0), vec3(1.0), i/4.0);
+	if(f < 4.0) return mix(vec3(1.0, 1.0, 0.0), vec3(1.0), i/4.0);
+	if(f < 5.0) return mix(vec3(1.0, 0.0, 1.0), vec3(1.0), i/4.0);
+	if(f < 6.0) return mix(vec3(0.0, 1.0, 1.0), vec3(1.0), i/4.0);
+	if(f < 7.0) return mix(vec3(1.0, 0.5, 0.0), vec3(1.0), i/4.0);
+	if(f < 8.0) return mix(vec3(1.0, 1.0, 1.0), vec3(1.0), i/4.0);
+	return vec3(1.0);
+}
+
+void main()
+{
+	ivec2 pixelCoord = ivec2(floor(texCoord.xy * vec2(textureSize(tex, 0).xy)));
+	float f = texelFetch(tex, pixelCoord, 0).x;
+	fragColor = vec4(color_map(f), 1.0);
+}

src/shaders/blit_div_vertex.glsl

@@ -0,0 +1,14 @@
+#version 300 es
+
+precision highp float;
+
+in vec2 pos;
+out vec2 texCoord;
+
+uniform mat4 mvp;
+
+void main()
+{
+	texCoord = 0.5*(pos + vec2(1,1));
+	gl_Position = mvp * vec4(pos, 0, 1);
+}

src/shaders/blit_residue_fragment.glsl

@@ -0,0 +1,65 @@
+#version 300 es
+
+precision highp float;
+precision highp sampler2D;
+
+in vec2 texCoord;
+out vec4 fragColor;
+
+uniform sampler2D xTex;
+uniform sampler2D bTex;
+
+float x(ivec2 coord)
+{
+	if(  coord.x <= 0
+	  || coord.x >= textureSize(xTex, 0).x
+	  || coord.y <= 0
+	  || coord.y >= textureSize(xTex, 0).y)
+	{
+		return(0.);
+	}
+	return(texelFetch(xTex, coord, 0).x);
+}
+
+float b(ivec2 coord)
+{
+	if(  coord.x <= 0
+	  || coord.x >= textureSize(bTex, 0).x
+	  || coord.y <= 0
+	  || coord.y >= textureSize(bTex, 0).y)
+	{
+		return(0.);
+	}
+	return(texelFetch(bTex, coord, 0).x);
+}
+
+vec3 color_map(float v)
+{
+	float logv = log(abs(v))/log(10.0);
+	float f = floor(logv + 7.0);
+	float i = floor(4.0*(logv + 7.0 - f));
+
+	if(f < 0.0) return vec3(0.0);
+	if(f < 1.0) return mix(vec3(1.0, 0.0, 0.0), vec3(1.0), i/4.0);
+	if(f < 2.0) return mix(vec3(0.0, 1.0, 0.0), vec3(1.0), i/4.0);
+	if(f < 3.0) return mix(vec3(0.0, 0.0, 1.0), vec3(1.0), i/4.0);
+	if(f < 4.0) return mix(vec3(1.0, 1.0, 0.0), vec3(1.0), i/4.0);
+	if(f < 5.0) return mix(vec3(1.0, 0.0, 1.0), vec3(1.0), i/4.0);
+	if(f < 6.0) return mix(vec3(0.0, 1.0, 1.0), vec3(1.0), i/4.0);
+	if(f < 7.0) return mix(vec3(1.0, 0.5, 0.0), vec3(1.0), i/4.0);
+	if(f < 8.0) return mix(vec3(1.0, 1.0, 1.0), vec3(1.0), i/4.0);
+	return vec3(1.0);
+}
+
+void main()
+{
+	ivec2 pixelCoord = ivec2(floor(texCoord.xy * vec2(textureSize(xTex, 0).xy)));
+
+	float tl = x(pixelCoord + ivec2(-1, 1));
+	float tr = x(pixelCoord + ivec2(1, 1));
+	float bl = x(pixelCoord + ivec2(-1, -1));
+	float br = x(pixelCoord + ivec2(1, -1));
+
+	float residue = b(pixelCoord) - (-tl - tr - bl - br + 4.*x(pixelCoord));
+	fragColor = vec4(color_map(residue), 1);
+}

src/shaders/blit_vertex.glsl

@@ -6,9 +6,13 @@ in vec2 pos;
 out vec2 texCoord;
 
 uniform mat4 mvp;
+uniform ivec2 gridSize;
 
 void main()
 {
-	texCoord = 0.5*(pos + vec2(1,1));
+	float margin = 32.;
+	float ratio = 1. - 2.*margin/float(gridSize.x);
+
+	texCoord = margin/float(gridSize.x) + ratio*(0.5*(pos + vec2(1,1)));
 	gl_Position = mvp * vec4(pos, 0, 1);
 }

src/shaders/divergence.glsl

@@ -8,36 +8,34 @@ out vec4 fragColor;
 
 uniform sampler2D src;
 
-float ux(ivec2 coord)
+vec2 u(ivec2 coord)
 {
-	return(texelFetch(src, coord, 0).x);
-}
-
-float uy(ivec2 coord)
-{
-	return(texelFetch(src, coord, 0).y);
+	return(texelFetch(src, coord, 0).xy);
 }
 
 void main()
 {
 	ivec2 pixelCoord = ivec2(floor(gl_FragCoord.xy));
 
-	ivec2 vr = pixelCoord + ivec2(1, 0);
-	ivec2 vl = pixelCoord - ivec2(1, 0);
-	ivec2 vt = pixelCoord + ivec2(0, 1);
-	ivec2 vb = pixelCoord - ivec2(0, 1);
-
-	float r = ux(vr);
-	float l = ux(vl);
-	float t = uy(vt);
-	float b = uy(vb);
-
-/*
-	vec2 c = texelFetch(src, pixelCoord, 0).xy;
-	if(vr.x >= textureSize(src, 0).x) { r = -c.x; }
-	if(vl.x < 0) { l = -c.x; }
-	if(vt.y >= textureSize(src, 0).y) { t = -c.y; }
-	if(vb.y < 0) { b = -c.y; }
-*/
-	fragColor = vec4(-(r - l + t - b)/2., 0, 0, 1);
+	if(  pixelCoord.x <= 0
+	  || pixelCoord.x >= textureSize(src, 0).x
+	  || pixelCoord.y <= 0
+	  || pixelCoord.y >= textureSize(src, 0).y)
+	{
+		fragColor = vec4(0, 0, 0, 1);
+	}
+	else
+	{
+		vec2 tl = u(pixelCoord + ivec2(-1, 0));
+		vec2 tr = u(pixelCoord);
+		vec2 bl = u(pixelCoord + ivec2(-1, -1));
+		vec2 br = u(pixelCoord + ivec2(0, -1));
+
+		float r = (tr.x + br.x)/2.;
+		float l = (tl.x + bl.x)/2.;
+		float t = (tl.y + tr.y)/2.;
+		float b = (bl.y + br.y)/2.;
+
+		fragColor = vec4(-2.*(r - l + t - b), 0, 0, 1);
+	}
 }

src/shaders/jacobi_step.glsl

@@ -12,9 +12,9 @@ uniform sampler2D bTex;
 float x(ivec2 coord)
 {
 	if(  coord.x <= 0
-	  || coord.x > textureSize(xTex, 0).x
+	  || coord.x >= textureSize(xTex, 0).x
 	  || coord.y <= 0
-	  || coord.y > textureSize(xTex, 0).y)
+	  || coord.y >= textureSize(xTex, 0).y)
 	{
 		return(0.);
 	}
@@ -23,6 +23,13 @@ float x(ivec2 coord)
 
 float b(ivec2 coord)
 {
+	if(  coord.x <= 0
+	  || coord.x >= textureSize(bTex, 0).x
+	  || coord.y <= 0
+	  || coord.y >= textureSize(bTex, 0).y)
+	{
+		return(0.);
+	}
 	return(texelFetch(bTex, coord, 0).x);
 }
 
@@ -30,14 +37,19 @@ void main()
 {
 	ivec2 pixelCoord = ivec2(floor(gl_FragCoord.xy));
 
-	ivec2 vr = pixelCoord + ivec2(1, 0);
-	ivec2 vl = pixelCoord - ivec2(1, 0);
-	ivec2 vt = pixelCoord + ivec2(0, 1);
-	ivec2 vb = pixelCoord - ivec2(0, 1);
-
-	float w = 0.8;
-	float standardJacobi = (x(vl) + x(vr) + x(vt) + x(vb) + b(pixelCoord))/4.;
-	float weightedJacobi = w*standardJacobi + (1.-w)*x(pixelCoord);
+	if(  pixelCoord.x <= 0
+	  || pixelCoord.y <= 0)
+	{
+		fragColor = vec4(0, 0, 0, 1);
+	}
+	else
+	{
+		float tl = x(pixelCoord + ivec2(-1, 1));
+		float tr = x(pixelCoord + ivec2(1, 1));
+		float bl = x(pixelCoord + ivec2(-1, -1));
+		float br = x(pixelCoord + ivec2(1, -1));
 
-	fragColor = vec4(weightedJacobi, 0, 0, 1);
+		float jacobi = (tl + tr + bl + br + b(pixelCoord))/4.;
+		fragColor = vec4(jacobi, 0, 0, 1);
+	}
 }

src/shaders/multigrid_correct.glsl

@@ -12,9 +12,28 @@ uniform float invGridSize;
 
 float e(ivec2 coord)
 {
+	if(  coord.x <= 0
+	  || coord.x >= textureSize(error, 0).x
+	  || coord.y <= 0
+	  || coord.y >= textureSize(error, 0).y)
+	{
+		return(0.);
+	}
 	return(texelFetch(error, coord, 0).x);
 }
 
+float p(ivec2 coord)
+{
+	if(  coord.x <= 0
+	  || coord.x >= textureSize(src, 0).x
+	  || coord.y <= 0
+	  || coord.y >= textureSize(src, 0).y)
+	{
+		return(0.);
+	}
+	return(texelFetch(src, coord, 0).x);
+}
+
 void main()
 {
 	ivec2 pixelCoord = ivec2(floor(gl_FragCoord.xy));
@@ -30,5 +49,5 @@ void main()
 	float bottomLerp = (1.-offset.x)*e(bl) + offset.x*e(br);
 	float bilerpError = (1.-offset.y)*bottomLerp + offset.y*topLerp;
 
-	fragColor = vec4(texture(src, texCoord).x + bilerpError, 0, 0, 1);
+	fragColor = vec4(p(pixelCoord) + bilerpError, 0, 0, 1);
 }

src/shaders/multigrid_restrict_residual.glsl

@@ -12,9 +12,9 @@ uniform sampler2D bTex;
 float x(ivec2 coord)
 {
 	if(  coord.x <= 0
-	  || coord.x > textureSize(xTex, 0).x
+	  || coord.x >= textureSize(xTex, 0).x
 	  || coord.y <= 0
-	  || coord.y > textureSize(xTex, 0).y)
+	  || coord.y >= textureSize(xTex, 0).y)
 	{
 		return(0.);
 	}
@@ -24,9 +24,9 @@ float x(ivec2 coord)
 float b(ivec2 coord)
 {
 	if(  coord.x <= 0
-	  || coord.x > textureSize(bTex, 0).x
+	  || coord.x >= textureSize(bTex, 0).x
 	  || coord.y <= 0
-	  || coord.y > textureSize(bTex, 0).y)
+	  || coord.y >= textureSize(bTex, 0).y)
 	{
 		return(0.);
 	}

src/shaders/splat.glsl

@@ -13,32 +13,16 @@ uniform float radius;
 uniform float additive;
 uniform float blending;
 
+uniform float randomize;
+
 void main()
 {
 	float d2 = dot(texCoord - splatPos, texCoord - splatPos);
 	float intensity = exp(-10.*d2/radius);
-
-	float n = 119.*gl_FragCoord.x + 107.*gl_FragCoord.y + 113.;
-	n = fract(n*fract(n/2.7182818));
-	intensity += intensity*abs(cos(2.*3.14159*n));
-
-	//NOTE: the force is added in the additive part.
 	vec2 force = splatColor.xy;
-	force.x += 0.2*sin(2.*3.14159*fract(n/1.144213));
-	force.y += 0.2*cos(2.*3.14159*fract(n/1.144213));
-
-	{
-		vec2 u = force/length(force);
-
-		vec2 v = texCoord - splatPos;
-		vec2 norm = v - dot(v, u)*u;
-		float dFromAxis = length(norm);
-
-		force += 30.*norm;
- 	}
 
 	vec3 u = texture(src, texCoord).xyz;
-	vec3 uAdd = u + intensity*vec3(force, 0);
+	vec3 uAdd = u + intensity*splatColor.xyz;
 	vec3 uBlend = u*(1.-intensity) + intensity * splatColor;
 
 	fragColor = vec4(uAdd*additive + uBlend*blending, 1);

src/shaders/subtract_pressure.glsl

@@ -17,11 +17,13 @@ vec2 u(ivec2 coord)
 
 float p(ivec2 coord)
 {
-	if(coord.x < 0) { coord.x = 0; }
-	if(coord.x >= textureSize(pressure, 0).x) { coord.x = textureSize(pressure, 0).x-1; }
-	if(coord.y < 0) { coord.y = 0; }
-	if(coord.y >= textureSize(pressure, 0).y) { coord.y = textureSize(pressure, 0).y-1; }
-
+	if(  coord.x <= 0
+	  || coord.x >= textureSize(pressure, 0).x
+	  || coord.y <= 0
+	  || coord.y >= textureSize(pressure, 0).y)
+	{
+		return(0.);
+	}
 	return(texelFetch(pressure, coord, 0).x);
 }
 
@@ -29,12 +31,17 @@ void main()
 {
 	ivec2 pixelCoord = ivec2(floor(gl_FragCoord.xy));
 
-	ivec2 vr = pixelCoord + ivec2(1, 0);
-	ivec2 vl = pixelCoord - ivec2(1, 0);
-	ivec2 vt = pixelCoord + ivec2(0, 1);
-	ivec2 vb = pixelCoord - ivec2(0, 1);
+	float tl = p(pixelCoord + ivec2(0, 1));
+	float tr = p(pixelCoord + ivec2(1, 1));
+	float bl = p(pixelCoord);
+	float br = p(pixelCoord + ivec2(1, 0));
+
+	float r = (tr + br)/2.;
+	float l = (tl + bl)/2.;
+	float t = (tl + tr)/2.;
+	float b = (bl + br)/2.;
 
-	vec2 gradP = vec2(p(vr) - p(vl), p(vt) - p(vb))/2.;
+	vec2 gradP = vec2(r - l, t - b);
 
 	fragColor = vec4(u(pixelCoord) - gradP, 0, 1);
 }