Spring-mass GLSL

This spring-mass system was adapted to the GPU in GLSL from Paul Bourke's example. I'm using it to sonify the fluid in a way that allows for more structure and flexibility than directly sampling the velocity field. I'll post more about this technique in the future.

The system uses a 1024 x 4 texture map with simple one-dimensional spring connections that bind each pixel to its horizontal neighbor, forming 4 lines. I extended Paul's example to include an anchor force which pulls the particles back to a rest position (in this case concentric rings), and an impulse force that applies velocity from an outside source.

layout(location=0) out vec3 fragColorPosition;
layout(location=1) out vec3 fragColorVelocity;

uniform float dt;
uniform float springConstant;
uniform float restLength;
uniform float dampingConstant;
uniform float viscousDrag;
uniform float anchorForce;
uniform float impulseForce;
uniform float mass;
uniform float fluidLengthThreshold;

vec2 springs[2] = vec2[2](vec2(-1.,0.), vec2(1.,0.));

const int POSITION_OLD = 0;
const int VELOCITY_OLD = 1;
const int VELOCITY_NEW = 2;
const int ANCHOR = 3;

const float smallf = 0.000001;

void main() {

    // Calculate forces
    vec3 pos0 = texture(sTD2DInputs[POSITION_OLD], vUV.st).xyz;
    vec3 vel0 = texture(sTD2DInputs[VELOCITY_OLD], vUV.st).xyz;
    vec3 vel_new = texture(sTD2DInputs[VELOCITY_NEW], pos0.xy).xyz; 
    vec3 anchor = texture(sTD2DInputs[ANCHOR], vUV.st).xyz;
    vec3 force = vec3(0.);

    force -= viscousDrag * vel0;
    force += vel_new * dt * dt * impulseForce;
    force += (anchor-pos0) * anchorForce;

    // spring interaction
    // 1d horizontal connection
    float cellSize = uTD2DInfos[POSITION_OLD].res.x; // 1/width

    for (int i = 0; i < 2; i++) {
        vec2 offset = vec2(cellSize, 0.);
        offset *= springs[i];
        vec2 coord = vUV.st + offset;
        
        if (coord.x < 0.) coord.x = 1.-cellSize;
        if (coord.x > 1.) coord.x = 0.+cellSize;

        vec3 pos1 = texture(sTD2DInputs[POSITION_OLD], coord).xyz;
        vec3 vel1 = texture(sTD2DInputs[VELOCITY_OLD], coord).xyz;

        vec3 dx = pos0 - pos1;
        float len = length(dx) + smallf;

        vec3 f = vec3(springConstant * (len - restLength)); 
        f += dampingConstant * (vel0 - vel1) * dx/vec3(len); 
        f *= -dx/vec3(len);

        force += f;
    }

    // Apply derivative
    vec3 dpdt = vel0;
    vec3 dvdt = force/vec3(mass);

    pos0 += dpdt * dt;
    vel0 += dvdt * dt;

    fragColorPosition = pos0;
    fragColorVelocity = vel0; 

} //main