+ if (num_nonzero_phases == 0) {
+ // Simplified version of the code below, assuming diff_flow == 0.0f everywhere.
+ diff_flow.x = omega * b.x * inv_A11;
+ diff_flow.y = omega * b.y * inv_A22;
+ } else {
+ // Subtract the missing terms from the right-hand side
+ // (it couldn't be done earlier, because we didn't know
+ // the values of the neighboring pixels; they change for
+ // each SOR iteration).
+ float smooth_l = zero_if_outside_border(texture(diffusivity_tex, tc_left));
+ float smooth_r = zero_if_outside_border(textureOffset(diffusivity_tex, tc_left, ivec2(1, 0)));
+ float smooth_d = zero_if_outside_border(texture(diffusivity_tex, tc_down));
+ float smooth_u = zero_if_outside_border(textureOffset(diffusivity_tex, tc_down, ivec2(0, 1)));
+ b += smooth_l * textureOffset(diff_flow_tex, tc, ivec2(-1, 0)).xy;
+ b += smooth_r * textureOffset(diff_flow_tex, tc, ivec2( 1, 0)).xy;
+ b += smooth_d * textureOffset(diff_flow_tex, tc, ivec2( 0, -1)).xy;
+ b += smooth_u * textureOffset(diff_flow_tex, tc, ivec2( 0, 1)).xy;
+
+ if (num_nonzero_phases == 1) {
+ diff_flow = vec2(0.0f);
+ } else {
+ diff_flow = texture(diff_flow_tex, tc).xy;
+ }
+
+ // From https://en.wikipedia.org/wiki/Successive_over-relaxation.
+ float sigma_u = A12 * diff_flow.y;
+ diff_flow.x += omega * ((b.x - sigma_u) * inv_A11 - diff_flow.x);
+ float sigma_v = A12 * diff_flow.x;
+ diff_flow.y += omega * ((b.y - sigma_v) * inv_A22 - diff_flow.y);
+ }