substep1.comp

#version 430 core
#extension GL_ARB_shading_language_include: require
#ifdef GL_NV_shader_atomic_float
#extension GL_NV_shader_atomic_float: enable
#endif
#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;
  float stress = -dt * 4 * E * p_vol * (J[p] - 1) / sqr(dx);
  mat2 affine = mat2(stress, 0, 0, stress) + p_mass * C[p];
  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 gvel = p_mass * v[p] + affine * dpos;
#ifdef GL_NV_shader_atomic_float
      atomicAdd(grid_v[g_].x, weight * gvel.x);
      atomicAdd(grid_v[g_].y, weight * gvel.y);
      atomicAdd(grid_m[g_], weight * p_mass);
#else
      float rhs;
      int old, new;
      rhs = weight * gvel.x;
      do {
        old = grid_v_i[g_].x;
        new = floatBitsToInt(intBitsToFloat(old) + rhs);
      } while (old != atomicCompSwap(grid_v_i[g_].x, old, new));
      rhs = weight * gvel.y;
      do {
        old = grid_v_i[g_].y;
        new = floatBitsToInt(intBitsToFloat(old) + rhs);
      } while (old != atomicCompSwap(grid_v_i[g_].y, old, new));
      rhs = weight * p_mass;
      do {
        old = grid_m_i[g_];
        new = floatBitsToInt(intBitsToFloat(old) + rhs);
      } while (old != atomicCompSwap(grid_m_i[g_], old, new));
#endif
    }
  }
}