orca/samples/fluid/src/main.c

/************************************************************//**
*
*	@file: main.cpp
*	@author: Martin Fouilleul
*	@date: 27/02/2022
*	@revision:
*
*****************************************************************/

#include"orca.h"
#include"glsl_shaders.h"

//----------------------------------------------------------------
//NOTE(martin): GL vertex struct and identifiers
//----------------------------------------------------------------
typedef struct Vertex { float x, y; } Vertex;

typedef struct advect_program
{
	GLuint prog;

	GLint pos;
	GLint src;
	GLint velocity;
	GLint delta;
	GLint dissipation;

} advect_program;

typedef struct div_program
{
	GLuint prog;
	GLint pos;
	GLint src;

} div_program;

typedef struct jacobi_program
{
	GLuint prog;
	GLint pos;
	GLint xTex;
	GLint bTex;

} 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;
	GLint pos;
	GLint xTex;
	GLint bTex;

} multigrid_restrict_residual_program;

typedef struct multigrid_correct_program
{
	GLuint prog;
	GLint pos;
	GLint src;
	GLint error;
	GLint invGridSize;

} multigrid_correct_program;

typedef struct subtract_program
{
	GLuint prog;

	GLint pos;
	GLint src;
	GLint pressure;
	GLint invGridSize;

} subtract_program;

typedef struct blit_program
{
	GLuint prog;

	GLint pos;
	GLint mvp;
	GLint gridSize;
	GLint tex;
} blit_program;

typedef struct splat_program
{
	GLuint prog;

	GLint pos;
	GLint src;
	GLint splatPos;
	GLint splatColor;
	GLint radius;
	GLint additive;
	GLint blending;
	GLint randomize;

} splat_program;

typedef struct frame_buffer
{
	GLuint textures[2];
	GLuint fbos[2];
} frame_buffer;

advect_program advectProgram;
div_program divProgram;
jacobi_program jacobiProgram;
multigrid_restrict_residual_program multigridRestrictResidualProgram;
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;

const int MULTIGRID_COUNT = 4;
frame_buffer pressureBuffer[4];
frame_buffer divBuffer[4];

GLuint vertexBuffer;

//----------------------------------------------------------------
//NOTE(martin): initialization
//----------------------------------------------------------------

GLuint compile_shader(const char* vs, const char* fs)
{
	GLuint vertexShader = glCreateShader(GL_VERTEX_SHADER);
	glShaderSource(vertexShader, 1, &vs, 0);
	glCompileShader(vertexShader);

	GLuint fragmentShader = glCreateShader(GL_FRAGMENT_SHADER);
	glShaderSource(fragmentShader, 1, &fs, 0);
	glCompileShader(fragmentShader);

	GLuint prog = glCreateProgram();
	glAttachShader(prog, vertexShader);
	glAttachShader(prog, fragmentShader);
	glLinkProgram(prog);


	//TODO errors
	int status = 0;
	glGetProgramiv(prog, GL_LINK_STATUS, &status);
	if(status != GL_TRUE)
	{
		log_error("program failed to link: ");
		int logSize = 0;
		glGetProgramiv(prog, GL_INFO_LOG_LENGTH, &logSize);

		mem_arena_scope scratch = mem_scratch_begin();
		char* log = mem_arena_alloc(scratch.arena, logSize);

		glGetProgramInfoLog(prog, logSize, 0, log);
		log_error("%s\n", log);

		mem_scratch_end(scratch);
	}

	int err = glGetError();
	if(err)
	{
		log_error("gl error %i\n", err);
	}

	return(prog);
}

void init_advect(advect_program* program)
{
	log_info("compiling advect...");
	program->prog = compile_shader(glsl_common_vertex, glsl_advect);
	program->pos = glGetAttribLocation(program->prog, "pos");
	program->src = glGetUniformLocation(program->prog, "src");
	program->velocity = glGetUniformLocation(program->prog, "velocity");
	program->delta = glGetUniformLocation(program->prog, "delta");
	program->dissipation = glGetUniformLocation(program->prog, "dissipation");
}

void init_div(div_program* program)
{
	log_info("compiling div...");
	program->prog = compile_shader(glsl_common_vertex, glsl_divergence);
	program->pos = glGetAttribLocation(program->prog, "pos");
	program->src = glGetUniformLocation(program->prog, "src");
}

void init_jacobi(jacobi_program* program)
{
	log_info("compiling jacobi...");
	program->prog = compile_shader(glsl_common_vertex, glsl_jacobi_step);
	program->pos = glGetAttribLocation(program->prog, "pos");
	program->xTex = glGetUniformLocation(program->prog, "xTex");
	program->bTex = glGetUniformLocation(program->prog, "bTex");
}

void init_multigrid_restrict_residual(multigrid_restrict_residual_program* program)
{
	log_info("compiling multigrid restrict residual...");
	program->prog = compile_shader(glsl_common_vertex, glsl_multigrid_restrict_residual);
	program->pos = glGetAttribLocation(program->prog, "pos");
	program->xTex = glGetUniformLocation(program->prog, "xTex");
	program->bTex = glGetUniformLocation(program->prog, "bTex");
}

void init_multigrid_correct(multigrid_correct_program* program)
{
	log_info("compiling multigrid correct...");
	program->prog = compile_shader(glsl_common_vertex, glsl_multigrid_correct);
	program->pos = glGetAttribLocation(program->prog, "pos");
	program->src = glGetUniformLocation(program->prog, "src");
	program->error = glGetUniformLocation(program->prog, "error");
	program->invGridSize = glGetUniformLocation(program->prog, "invGridSize");
}

void init_subtract(subtract_program* program)
{
	log_info("compiling subtract...");
	program->prog = compile_shader(glsl_common_vertex, glsl_subtract_pressure);
	program->pos = glGetAttribLocation(program->prog, "pos");
	program->src = glGetUniformLocation(program->prog, "src");
	program->pressure = glGetUniformLocation(program->prog, "pressure");
	program->invGridSize = glGetUniformLocation(program->prog, "invGridSize");
}

void init_splat(splat_program* program)
{
	log_info("compiling splat...");
	program->prog = compile_shader(glsl_common_vertex, glsl_splat);
	program->pos = glGetAttribLocation(program->prog, "pos");
	program->src = glGetUniformLocation(program->prog, "src");
	program->splatPos = glGetUniformLocation(program->prog, "splatPos");
	program->splatColor = glGetUniformLocation(program->prog, "splatColor");
	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)
{
	log_info("compiling blit...");
	program->prog = compile_shader(glsl_blit_vertex, glsl_blit_fragment);
	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)
{
	log_info("compiling blit div...");
	program->prog = compile_shader(glsl_blit_div_vertex, glsl_blit_div_fragment);
	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)
{
	log_info("compiling blit residue...");
	program->prog = compile_shader(glsl_blit_div_vertex, glsl_blit_residue_fragment);
	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;
	glGenTextures(1, &texture);
	glBindTexture(GL_TEXTURE_2D, texture);
	glTexImage2D(GL_TEXTURE_2D, 0, internalFormat, width, height, 0, format, type, initData);
	glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
	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);
}

GLuint create_fbo(GLuint texture)
{
	GLuint fbo;
	glGenFramebuffers(1, &fbo);
	glBindFramebuffer(GL_FRAMEBUFFER, fbo);
	glBindTexture(GL_TEXTURE_2D, texture);
	glFramebufferTexture2D(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT0, GL_TEXTURE_2D, texture, 0);
	return(fbo);
}

void init_frame_buffer(frame_buffer* framebuffer,
				       int width,
                       int height,
                       GLenum internalFormat,
                       GLenum format,
                       GLenum type,
                       char* initData)
{
	for(int i=0; i<2; i++)
	{
		framebuffer->textures[i] = create_texture(width, height, internalFormat, format, type, initData);
		framebuffer->fbos[i] = create_fbo(framebuffer->textures[i]);
	}

	GLenum err = glCheckFramebufferStatus(GL_FRAMEBUFFER);
	if(err != GL_FRAMEBUFFER_COMPLETE)
	{
		log_info("Frame buffer incomplete, %i", err);
	}
}

void frame_buffer_swap(frame_buffer* buffer)
{
	GLuint tmp = buffer->fbos[0];
	buffer->fbos[0] = buffer->fbos[1];
	buffer->fbos[1] = tmp;

	tmp = buffer->textures[0];
	buffer->textures[0] = buffer->textures[1];
	buffer->textures[1] = tmp;
}

//----------------------------------------------------------------
//NOTE(martin): entry point
//----------------------------------------------------------------

#define texWidth (256)
#define texHeight (256)

float colorInitData[texWidth][texHeight][4] = {0};
float velocityInitData[texWidth][texHeight][4] = {0};

const float EPSILON = 1.,
            INV_GRID_SIZE = 1./(float)texWidth,
            DELTA = 1./120.;

const GLenum TEX_INTERNAL_FORMAT = GL_RGBA32F;
const GLenum TEX_FORMAT = GL_RGBA;
const GLenum TEX_TYPE = GL_FLOAT;

#define square(x) ((x)*(x))

/*
void reset_texture(GLuint texture, float width, float height, char* initData)
{
	glBindTexture(GL_TEXTURE_2D, texture);
	glTexImage2D(GL_TEXTURE_2D, 0, TEX_INTERNAL_FORMAT, width, height, 0, TEX_FORMAT, TEX_TYPE, initData);
}

static bool resetCmd = false;

void reset()
{
//	resetCmd = true;
	log_info("reset");

	reset_texture(colorBuffer.textures[0], texWidth, texHeight, (char*)colorInitData);
	reset_texture(colorBuffer.textures[1], texWidth, texHeight, (char*)colorInitData);
	reset_texture(velocityBuffer.textures[0], texWidth, texHeight, (char*)velocityInitData);
	reset_texture(velocityBuffer.textures[1], texWidth, texHeight, (char*)velocityInitData);

	int gridFactor = 1;
	for(int i=0; i<MULTIGRID_COUNT; i++)
	{
		reset_texture(pressureBuffer[i].textures[0], texWidth/gridFactor, texHeight/gridFactor, 0);
		reset_texture(pressureBuffer[i].textures[1], texWidth/gridFactor, texHeight/gridFactor, 0);

		gridFactor *= 2;
	}
}

*/

typedef struct mouse_input
{
	float x;
	float y;
	float deltaX;
	float deltaY;
	bool down;

} mouse_input;

mouse_input mouseInput = {0};

int frameWidth = 800;
int frameHeight = 600;


ORCA_EXPORT void OnMouseDown(int button)
{
	mouseInput.down = true;
}

ORCA_EXPORT void OnMouseUp(int button)
{
	mouseInput.down = false;
}

ORCA_EXPORT void OnMouseMove(float x, float y, float dx, float dy)
{
	mouseInput.x = x * 2;
	mouseInput.y = y * 2;
	mouseInput.deltaX = dx * 2;
	mouseInput.deltaY = dy * 2;
}

void init_color_checker()
{
	for(int i=0; i<texHeight; i++)
	{
		for(int j=0; j<texWidth; j++)
		{
			float u = j/(float)texWidth;
			float v = i/(float)texWidth;
			float value = ((int)(u*10)%2) == ((int)(v*10)%2) ? 1. : 0.;

			for(int k = 0; k<3; k++)
			{
				colorInitData[i][j][k] = value;
			}
			colorInitData[i][j][3] = 1;
		}
	}
}

void init_velocity_vortex()
{
	for(int i=0; i<texHeight; i++)
	{
		for(int j=0; j<texWidth; j++)
		{
			float x = 2*j/(float)texWidth - 1;
			float y = 2*i/(float)texWidth - 1;
			velocityInitData[i][j][0] = sinf(2*M_PI*y);
			velocityInitData[i][j][1] = sinf(2*M_PI*x);

		}
	}
}

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]);

	glActiveTexture(GL_TEXTURE0);
	glBindTexture(GL_TEXTURE_2D, velocityBuffer.textures[0]);
	glUniform1i(splatProgram.src, 0);

	glUniform2f(splatProgram.splatPos, splatPosX, splatPosY);
	glUniform3f(splatProgram.splatColor, splatVelX, splatVelY, 0);
	glUniform1f(splatProgram.additive, 1);
	glUniform1f(splatProgram.blending, 0);

	glUniform1f(splatProgram.radius, radius);

	glDrawArrays(GL_TRIANGLES, 0, 6);

	frame_buffer_swap(&velocityBuffer);

	// dye
	glBindFramebuffer(GL_FRAMEBUFFER, colorBuffer.fbos[1]);

	glActiveTexture(GL_TEXTURE0);
	glBindTexture(GL_TEXTURE_2D, colorBuffer.textures[0]);
	glUniform1i(splatProgram.src, 0);

	glUniform2f(splatProgram.splatPos, splatPosX, splatPosY);
	glUniform3f(splatProgram.splatColor, r, g, b);
	glUniform1f(splatProgram.additive, 0);
	glUniform1f(splatProgram.blending, 1);
	glUniform1f(splatProgram.radius, radius);

	glDrawArrays(GL_TRIANGLES, 0, 6);

	frame_buffer_swap(&colorBuffer);

}

void jacobi_solve(frame_buffer* x, frame_buffer* b, float invGridSize, int iterationCount)
{
	glUseProgram(jacobiProgram.prog);

	for(int i=0; i<iterationCount; i++)
	{
		glBindFramebuffer(GL_FRAMEBUFFER, x->fbos[1]);

		glActiveTexture(GL_TEXTURE0);
		glBindTexture(GL_TEXTURE_2D, x->textures[0]);
		glUniform1i(jacobiProgram.xTex, 0);

		glActiveTexture(GL_TEXTURE1);
		glBindTexture(GL_TEXTURE_2D, b->textures[0]);
		glUniform1i(jacobiProgram.bTex, 1);

		glDrawArrays(GL_TRIANGLES, 0, 6);

		frame_buffer_swap(x);
	}
}

void multigrid_coarsen_residual(frame_buffer* output, frame_buffer* x, frame_buffer* b, float invFineGridSize)
{
	//NOTE: compute residual and downsample to coarser grid, put result in coarser buffer
	glUseProgram(multigridRestrictResidualProgram.prog);
	glBindFramebuffer(GL_FRAMEBUFFER, output->fbos[1]);

	glActiveTexture(GL_TEXTURE0);
	glBindTexture(GL_TEXTURE_2D, x->textures[0]);
	glUniform1i(multigridRestrictResidualProgram.xTex, 0);

	glActiveTexture(GL_TEXTURE1);
	glBindTexture(GL_TEXTURE_2D, b->textures[0]);
	glUniform1i(multigridRestrictResidualProgram.bTex, 1);

	glDrawArrays(GL_TRIANGLES, 0, 6);

	frame_buffer_swap(output);
}

void multigrid_prolongate_and_correct(frame_buffer* x, frame_buffer* error, float invFineGridSize)
{
	//NOTE: correct finer pressure
	glUseProgram(multigridCorrectProgram.prog);
	glBindFramebuffer(GL_FRAMEBUFFER, x->fbos[1]);

	glActiveTexture(GL_TEXTURE0);
	glBindTexture(GL_TEXTURE_2D, x->textures[0]);
	glUniform1i(multigridCorrectProgram.src, 0);

	glActiveTexture(GL_TEXTURE1);
	glBindTexture(GL_TEXTURE_2D, error->textures[0]);
	glUniform1i(multigridCorrectProgram.error, 1);

	glUniform1f(multigridCorrectProgram.invGridSize, invFineGridSize);

	glDrawArrays(GL_TRIANGLES, 0, 6);

	frame_buffer_swap(x);
}

void multigrid_clear(frame_buffer* error)
{
	glBindFramebuffer(GL_FRAMEBUFFER, error->fbos[0]);
	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 offset = margin/texWidth;
		float ratio = 1 - 2*margin/texWidth;

		float splatPosX = (mouseInput.x/frameWidth)*ratio + offset;
		float splatPosY = (1 - mouseInput.y/frameHeight)*ratio + offset;

		float splatVelX = (10000.*DELTA*mouseInput.deltaX/frameWidth)*ratio;
		float splatVelY = (-10000.*DELTA*mouseInput.deltaY/frameWidth)*ratio;

		float intensity = 100*sqrtf(square(ratio*mouseInput.deltaX/frameWidth) + square(ratio*mouseInput.deltaY/frameHeight));

		float r = intensity * (sinf(2*M_PI*0.1*t) + 1);
		float g = 0.5*intensity * (cosf(2*M_PI*0.1/M_E*t + 654) + 1);
		float b = intensity * (sinf(2*M_PI*0.1/M_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;
	}
}

float testDiv[texWidth/2][texWidth/2][4];

mg_surface surface;

ORCA_EXPORT void OnInit()
{
	log_info("Hello, world (from C)");

	surface = mg_surface_gles();
	mg_surface_prepare(surface);

//	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
	log_info("create color buffer");
	init_frame_buffer(&colorBuffer, texWidth, texHeight, TEX_INTERNAL_FORMAT, TEX_FORMAT, TEX_TYPE, (char*)colorInitData);
	log_info("create velocity buffer");
	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++)
	{
		log_info("create div buffer %i", i);
		init_frame_buffer(&divBuffer[i], texWidth/gridFactor, texHeight/gridFactor, TEX_INTERNAL_FORMAT, TEX_FORMAT, TEX_TYPE, 0);
		log_info("create pressure buffer %i", i);
		init_frame_buffer(&pressureBuffer[i], texWidth/gridFactor, texHeight/gridFactor, TEX_INTERNAL_FORMAT, TEX_FORMAT, TEX_TYPE, 0);
		gridFactor *= 2;
	}

	// init vertex buffer
	static Vertex vertices[6] =
	{
		{-1, -1},
		{ 1, -1},
		{ 1,  1},
		{-1, -1},
		{ 1,  1},
		{ -1, 1}
	};

	//WARN: we assume blitProgram.pos == advectProgram.pos, is there a situation where it wouldn't be true??
	GLuint vertexBuffer = 0;
	glGenBuffers(1, &vertexBuffer);
	glBindBuffer(GL_ARRAY_BUFFER, vertexBuffer);
	glBufferData(GL_ARRAY_BUFFER, 6*sizeof(Vertex), vertices, GL_STATIC_DRAW);

    glEnableVertexAttribArray(0);
    glVertexAttribPointer(blitProgram.pos, 2, GL_FLOAT, GL_FALSE, 0, 0);


	for(int i=0; i<texWidth/2; i++)
	{
		for(int j=0; j<texHeight/2; j++)
		{
			testDiv[i][j][0] = 0.5 + 0.5*cosf(j/100.*3.14159 + i/100.*1.2139);
		}
	}
}

ORCA_EXPORT void OnFrameResize(u32 width, u32 height)
{
	frameWidth = width*2;
	frameHeight = height*2;
}

ORCA_EXPORT void OnFrameRefresh()
{
	float aspectRatio = texWidth/texHeight; //TODO replace with actual aspect ratio?

	static float t = 0;
	t += 1./60.;

	mg_surface_prepare(surface);

	glViewport(0, 0, texWidth, texHeight);

	//NOTE: advect velocity thru itself
	glUseProgram(advectProgram.prog);
	glBindFramebuffer(GL_FRAMEBUFFER, velocityBuffer.fbos[1]);

	glActiveTexture(GL_TEXTURE0);
	glBindTexture(GL_TEXTURE_2D, velocityBuffer.textures[0]);
	glUniform1i(advectProgram.src, 0);

	glActiveTexture(GL_TEXTURE1);
	glBindTexture(GL_TEXTURE_2D, velocityBuffer.textures[0]);
	glUniform1i(advectProgram.velocity, 1);

	glUniform1f(advectProgram.delta, DELTA);
	glUniform1f(advectProgram.dissipation, 0.01);

	glDrawArrays(GL_TRIANGLES, 0, 6);

	frame_buffer_swap(&velocityBuffer);

	/*
	//DEBUG
	static bool splatTrig = false;
	static bool splat = false;
	static float splatStart = 0;
	static int splatDir = 0;

	static int frameCount = 0;

	if(resetCmd)
	{
		frameCount = 0;
		splat = true;
		splatStart = frameT;
	}

	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;

	if(frameCount>20)
	{
		return;
	}
	frameCount++;
*/

	input_splat(t);


	//NOTE: compute divergence of advected velocity
	glUseProgram(divProgram.prog);
	glBindFramebuffer(GL_FRAMEBUFFER, divBuffer[0].fbos[1]);

	glActiveTexture(GL_TEXTURE0);
	glBindTexture(GL_TEXTURE_2D, velocityBuffer.textures[0]);
	glUniform1i(divProgram.src, 0);

	glDrawArrays(GL_TRIANGLES, 0, 6);

	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]);
		jacobi_solve(&pressureBuffer[0], &divBuffer[0], INV_GRID_SIZE, texWidth*texHeight);
	#else
		multigrid_clear(&pressureBuffer[0]);

		for(int i=0; i<1; i++)
		{
			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]);
			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]);
			jacobi_solve(&pressureBuffer[2], &divBuffer[2], 4*INV_GRID_SIZE, 30);

			multigrid_prolongate_and_correct(&pressureBuffer[1], &pressureBuffer[2], 2*INV_GRID_SIZE);
			jacobi_solve(&pressureBuffer[1], &divBuffer[1], 2*INV_GRID_SIZE, 8);

			multigrid_prolongate_and_correct(&pressureBuffer[0], &pressureBuffer[1], INV_GRID_SIZE);
			jacobi_solve(&pressureBuffer[0], &divBuffer[0], INV_GRID_SIZE, 4);
		}
	#endif




	//NOTE: subtract pressure gradient to advected velocity
	glUseProgram(subtractProgram.prog);
	glBindFramebuffer(GL_FRAMEBUFFER, velocityBuffer.fbos[1]);

	glActiveTexture(GL_TEXTURE0);
	glBindTexture(GL_TEXTURE_2D, velocityBuffer.textures[0]);
	glUniform1i(subtractProgram.src, 0);

	glActiveTexture(GL_TEXTURE1);
	glBindTexture(GL_TEXTURE_2D, pressureBuffer[0].textures[0]);
	glUniform1i(subtractProgram.pressure, 1);

	glUniform1f(subtractProgram.invGridSize, INV_GRID_SIZE);

	glDrawArrays(GL_TRIANGLES, 0, 6);

	frame_buffer_swap(&velocityBuffer);

	//NOTE: Advect color through corrected velocity field
	glUseProgram(advectProgram.prog);
	glBindFramebuffer(GL_FRAMEBUFFER, colorBuffer.fbos[1]);

	glActiveTexture(GL_TEXTURE0);
	glBindTexture(GL_TEXTURE_2D, colorBuffer.textures[0]);
	glUniform1i(advectProgram.src, 0);

	glActiveTexture(GL_TEXTURE1);
	glBindTexture(GL_TEXTURE_2D, velocityBuffer.textures[0]);
	glUniform1i(advectProgram.velocity, 1);

	glUniform1f(advectProgram.delta, DELTA);

	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, frameWidth, frameHeight);

	float displayMatrix[16] = {
		1/aspectRatio, 0, 0, 0,
		0, 1, 0, 0,
		0, 0, 1, 0,
		0, 0, 0, 1 };

/*
	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);

	glActiveTexture(GL_TEXTURE0);
	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: 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(divProgram.src, 0);

	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);

//*/

	mg_surface_present(surface);
}