substep3.comp

#version 430 core
#extension GL_ARB_shading_language_include: require
#include </glcommon.h>

layout (local_size_x = N1) in;

float sqr(float x) {
  return x * x;
}

float wfunc(int i, float fx) {
  if (i == 0) {
    return 0.5 * sqr(1.5 - fx);
  } else if (i == 1) {
    return 0.75 - sqr(fx - 1.0);
  } else {
    return 0.5 * sqr(fx - 0.5);
  }
}

void main()
{
  uint p = gl_GlobalInvocationID.x;
  vec2 Xp = x[p] / dx;
  ivec2 base = ivec2(Xp - 0.5);
  vec2 fx = Xp - base;
  vec2 new_v = vec2(0);
  mat2 new_C = mat2(0);
  for (int i = 0; i < 3; i++) {
    for (int j = 0; j < 3; j++) {
      ivec2 off = ivec2(i, j);
      vec2 dpos = (off - fx) * dx;
      float weight = wfunc(i, fx.x) * wfunc(j, fx.y);
      ivec2 gij = base + off;
      int g_ = gij.x * G + gij.y;
      vec2 g_v = grid_v[g_];
      new_v += weight * g_v;
      new_C += 4 * weight * outerProduct(g_v, dpos) / sqr(dx);
    }
  }
  v[p] = new_v;
  x[p] += dt * v[p];
  J[p] *= 1 + dt * (new_C[0][0] + new_C[1][1]);
  C[p] = new_C;
}