From 2767a75a1bd4450fd47c060eff18367e06c1532a Mon Sep 17 00:00:00 2001 From: Martin Fouilleul Date: Thu, 14 Jul 2022 18:28:50 +0200 Subject: [PATCH] vertex multigrid --- src/main.c | 271 +++++++++++++++++++++++-- src/shaders/advect.glsl | 12 +- src/shaders/blit_residue_fragment.glsl | 65 ++++++ src/shaders/divergence.glsl | 41 +--- src/shaders/jacobi_step.glsl | 19 +- src/shaders/splat.glsl | 23 +-- src/shaders/subtract_pressure.glsl | 34 +--- 7 files changed, 355 insertions(+), 110 deletions(-) create mode 100644 src/shaders/blit_residue_fragment.glsl diff --git a/src/main.c b/src/main.c index a2bfb4c..38c1bf8 100644 --- a/src/main.c +++ b/src/main.c @@ -25,6 +25,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) { @@ -134,6 +135,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; @@ -204,6 +215,7 @@ subtract_program subtractProgram; splat_program splatProgram; blit_program blitProgram; blit_program blitDivProgram; +blit_residue_program blitResidueProgram; frame_buffer colorBuffer; frame_buffer velocityBuffer; @@ -325,6 +337,16 @@ void init_blit_div(blit_program* program) 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_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) { @@ -384,12 +406,9 @@ void frame_buffer_swap(frame_buffer* buffer) //---------------------------------------------------------------- //NOTE(martin): entry point //---------------------------------------------------------------- -/* -#define texWidth 512 -#define texHeight 512 -*/ -#define texWidth (128) -#define texHeight (128) + +#define texWidth (256) +#define texHeight (256) float colorInitData[texWidth][texHeight][4] = {0}; float velocityInitData[texWidth][texHeight][4] = {0}; @@ -408,11 +427,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; +//static bool resetCmd = false; void reset() { - resetCmd = true; +// resetCmd = true; console_log_fmt("reset"); reset_texture(colorBuffer.textures[0], texWidth, texHeight, (char*)colorInitData); @@ -490,6 +509,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); + } } } @@ -504,8 +524,8 @@ int init(float canvasSize) WAJS_SetupCanvas(canvasSize, canvasSize); - init_color_checker(); - init_velocity_vortex(); +// init_color_checker(); +// init_velocity_vortex(); // init programs init_advect(&advectProgram); @@ -513,6 +533,7 @@ int init(float canvasSize) init_jacobi(&jacobiProgram); init_multigrid_restrict_residual(&multigridRestrictResidualProgram); init_multigrid_correct(&multigridCorrectProgram); + init_blit_residue(&blitResidueProgram); init_subtract(&subtractProgram); init_splat(&splatProgram); @@ -565,6 +586,56 @@ int init(float canvasSize) return(0); } +void apply_splat(float splatPosX, float splatPosY, 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, 0.01); + + 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, 0.01); + + 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); @@ -587,6 +658,79 @@ void jacobi_solve(frame_buffer* x, frame_buffer* b, float invGridSize, int itera } } +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)) + { + 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)); + + 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); + + apply_splat(splatPosX, splatPosY, splatVelX, splatVelY, r, g, b, false); + + mouseInput.deltaX = 0; + mouseInput.deltaY = 0; + } +} + // This function is called by loader.js every frame void WAFNDraw() { @@ -617,6 +761,61 @@ void WAFNDraw() 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]); @@ -633,7 +832,34 @@ void WAFNDraw() glBindFramebuffer(GL_FRAMEBUFFER, pressureBuffer[0].fbos[1]); glClear(GL_COLOR_BUFFER_BIT); - jacobi_solve(&pressureBuffer[0], &divBuffer[0], INV_GRID_SIZE, 20000); + #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<2; 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, 40); + + 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); @@ -667,7 +893,7 @@ void WAFNDraw() glUniform1f(advectProgram.delta, DELTA); - glUniform1f(advectProgram.dissipation, 0.3); + glUniform1f(advectProgram.dissipation, 0.001); glDrawArrays(GL_TRIANGLES, 0, 6); @@ -675,6 +901,9 @@ void WAFNDraw() //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, @@ -682,8 +911,22 @@ void WAFNDraw() 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); @@ -695,6 +938,7 @@ void WAFNDraw() glDrawArrays(GL_TRIANGLES, 0, 6); /*/ + //NOTE: recompute divergence of (corrected) velocity glUseProgram(divProgram.prog); glBindFramebuffer(GL_FRAMEBUFFER, divBuffer[0].fbos[1]); @@ -719,5 +963,6 @@ void WAFNDraw() glUniformMatrix4fv(blitDivProgram.mvp, 1, GL_FALSE, displayMatrix); glDrawArrays(GL_TRIANGLES, 0, 6); + //*/ } diff --git a/src/shaders/advect.glsl b/src/shaders/advect.glsl index 283a5ad..28cbd2c 100644 --- a/src/shaders/advect.glsl +++ b/src/shaders/advect.glsl @@ -13,22 +13,14 @@ uniform float dissipation; vec2 u(ivec2 coord) { - if( coord.x <= 0 - || coord.x >= textureSize(velocity, 0).x - || coord.y <= 0 - || coord.y >= textureSize(velocity, 0).y) - { - return(vec2(0.)); - } - return(texelFetch(velocity, coord, 0).xy); } vec4 q(ivec2 coord) { - if( coord.x <= 0 + if( coord.x < 0 || coord.x >= textureSize(src, 0).x - || coord.y <= 0 + || coord.y < 0 || coord.y >= textureSize(src, 0).y) { return(vec4(0.)); diff --git a/src/shaders/blit_residue_fragment.glsl b/src/shaders/blit_residue_fragment.glsl new file mode 100644 index 0000000..f630f41 --- /dev/null +++ b/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); +} diff --git a/src/shaders/divergence.glsl b/src/shaders/divergence.glsl index bf085bd..8b23781 100644 --- a/src/shaders/divergence.glsl +++ b/src/shaders/divergence.glsl @@ -8,28 +8,9 @@ out vec4 fragColor; uniform sampler2D src; -float ux(ivec2 coord) +vec2 u(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); -} - -float uy(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).y); + return(texelFetch(src, coord, 0).xy); } void main() @@ -45,16 +26,16 @@ void main() } else { - ivec2 vr = pixelCoord + ivec2(1, 0); - ivec2 vl = pixelCoord - ivec2(1, 0); - ivec2 vt = pixelCoord + ivec2(0, 1); - ivec2 vb = pixelCoord - ivec2(0, 1); + 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 = ux(vr); - float l = ux(vl); - float t = uy(vt); - float b = uy(vb); + 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(-(r - l + t - b)/2., 0, 0, 1); + fragColor = vec4(-2.*(r - l + t - b), 0, 0, 1); } } diff --git a/src/shaders/jacobi_step.glsl b/src/shaders/jacobi_step.glsl index 6493091..35418b9 100644 --- a/src/shaders/jacobi_step.glsl +++ b/src/shaders/jacobi_step.glsl @@ -38,23 +38,18 @@ void main() ivec2 pixelCoord = ivec2(floor(gl_FragCoord.xy)); if( pixelCoord.x <= 0 - || pixelCoord.x >= textureSize(xTex, 0).x - || pixelCoord.y <= 0 - || pixelCoord.y >= textureSize(xTex, 0).y) + || pixelCoord.y <= 0) { fragColor = vec4(0, 0, 0, 1); } else { - 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 = 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 w = 0.8; - float standardJacobi = (x(vl) + x(vr) + x(vt) + x(vb) + b(pixelCoord))/4.; - float weightedJacobi = w*standardJacobi + (1.-w)*x(pixelCoord); - - fragColor = vec4(weightedJacobi, 0, 0, 1); + float jacobi = (tl + tr + bl + br + b(pixelCoord))/4.; + fragColor = vec4(jacobi, 0, 0, 1); } } diff --git a/src/shaders/splat.glsl b/src/shaders/splat.glsl index ae76b3c..34a1347 100644 --- a/src/shaders/splat.glsl +++ b/src/shaders/splat.glsl @@ -21,29 +21,8 @@ void main() float intensity = exp(-10.*d2/radius); vec2 force = splatColor.xy; - if(randomize != 0.0) - { - 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. - - 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); diff --git a/src/shaders/subtract_pressure.glsl b/src/shaders/subtract_pressure.glsl index 7e6827d..7905df9 100644 --- a/src/shaders/subtract_pressure.glsl +++ b/src/shaders/subtract_pressure.glsl @@ -12,13 +12,6 @@ uniform float invGridSize; vec2 u(ivec2 coord) { - if( coord.x <= 0 - || coord.x >= textureSize(src, 0).x - || coord.y <= 0 - || coord.y >= textureSize(src, 0).y) - { - return(vec2(0.)); - } return(texelFetch(src, coord, 0).xy); } @@ -38,22 +31,17 @@ void main() { ivec2 pixelCoord = ivec2(floor(gl_FragCoord.xy)); - 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 - { - 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)); - vec2 gradP = vec2(p(vr) - p(vl), p(vt) - p(vb))/2.; + float r = (tr + br)/2.; + float l = (tl + bl)/2.; + float t = (tl + tr)/2.; + float b = (bl + br)/2.; - fragColor = vec4(u(pixelCoord) - gradP, 0, 1); - } + vec2 gradP = vec2(r - l, t - b); + + fragColor = vec4(u(pixelCoord) - gradP, 0, 1); }