initialize each vertex's normal vector to (0,0,0)
parallel for each face in mesh:
  for each vertex in face:
    atomic-add face's normal to vertex's normal // this is scatter
parallel for each vertex in mesh:
  normalize(vertex normal)

@group(0) @binding(0) var<storage> inputs: array<f32>;
@group(0) @binding(1) var<storage, read_write> output: f32;

// The cell holds u32 values.
var<workgroup> sum: atomic<u32>;

alias cell = ptr<workgroup,atomic<u32>>;

// Adds f32 'value' to the value in 'cell', atomically.
// Perform atomic compare-exchange with u32 type, and bitcast in and out of f32.
fn atomic_add(sum_cell: cell, value: f32) -> f32 {
  // Initializing to 0 forces second iteration in almost all cases.
  var old = 0u; //  alternately, atomicLoad(sum_cell);
  loop {
    let new_value = value + bitcast<f32>(old);
    let exchange_result = atomicCompareExchangeWeak(sum_cell, old, bitcast<u32>(new_value));
    if exchange_result.exchanged {
       return new_value;
    }
    old = exchange_result.old_value;
  }
}

@compute @workgroup_size(32)
fn main(@builtin(global_invocation_id) gid: vec3u) {
  atomic_add(&sum, inputs[gid.x]);
  workgroupBarrier();
  if gid.x == 0 {
    output = bitcast<f32>(atomicLoad(&sum));
  }
}

parallel for each vertex in mesh:
  initialize each vertex's normal vector to (0,0,0)
  for each face in all faces:
    if vertex in face:
      add face.normal to my normal
  normalize(vertex normal)

      const facet_normals_module = device.createShaderModule({
        label: "compute facet normals module",
        code: /* wgsl */ `
                    /* output */
                    @group(0) @binding(0) var<storage, read_write> facet_normals: array<vec3f>;
                    /* input */
                    @group(0) @binding(1) var<storage, read> vertices: array<vec3f>;
                    @group(0) @binding(2) var<storage, read> triangle_indices: array<u32>;

                     /** Algorithm:
                      * For tri in all triangles:
                      *   Fetch all 3 vertices of tri
                      *   Compute normalize(cross(v1-v0, v2-v0))
                      *   For each vertex in tri:
                      *     Atomically add it to vertex_normals[vertex]
                      *     /* Can't do this! No f32 atomics */
                      * For vertex in all vertices:
                      *   Normalize vertex_normals[vertex]
                      *
                      * OK, so we can't do this approach w/o f32 atomics
                      * So we will instead convert this scatter to gather
                      * This is wasteful; every vertex will walk the entire
                      *   index array looking for matches.
                      * Could alternately build a mapping of {vtx->facet}
                      *
                      * (1) For tri in all triangles:
                      *   Fetch all 3 vertices of tri
                      *   Compute normalize(cross(v1-v0, v2-v0))
                      *   Store that vector as a facet normal
                      * (2) For vertex in all vertices:
                      *   normal[vertex] = (0,0,0)
                      *   For tri in all triangles:
                      *     // note expensive doubly-nested loop!
                      *     if my vertex is in that triangle:
                      *       normal[vertex] += facet_normal[tri]
                      *   normalize(normal[vertex])
                      */
                    @compute @workgroup_size(${WORKGROUP_SIZE}) fn facet_normals_kernel(
                      @builtin(global_invocation_id) id: vec3u) {
                        let tri = id.x;
                        if (tri < arrayLength(&facet_normals)) {
                          /* note triangle_indices is u32 not vec3, do math accordingly */
                          let v0: vec3f = vertices[triangle_indices[tri * 3]];
                          let v1: vec3f = vertices[triangle_indices[tri * 3 + 1]];
                          let v2: vec3f = vertices[triangle_indices[tri * 3 + 2]];
                          facet_normals[tri] = normalize(cross(v1-v0, v2-v0));
                        }
                      }
                  `,
      });

      const vertex_normals_module = device.createShaderModule({
        label: "compute vertex normals module",
        code: /* wgsl */ `
                    /* output */
                    @group(0) @binding(0) var<storage, read_write> vertex_normals: array<vec3f>;
                    /* input */
                    @group(0) @binding(1) var<storage, read> facet_normals: array<vec3f>;
                    @group(0) @binding(2) var<storage, read> triangle_indices: array<u32>;

                    /* see facet_normals_module for algorithm */

                    @compute @workgroup_size(${WORKGROUP_SIZE}) fn vertex_normals_kernel(
                      @builtin(global_invocation_id) id: vec3u) {
                        let vtx = id.x;
                        if (vtx < arrayLength(&vertex_normals)) {
                          vertex_normals[vtx] = vec3f(0, 0, 0);
                          /* note triangle_indices is u32 not vec3, do math accordingly */
                          for (var tri: u32 = 0; tri < arrayLength(&triangle_indices) / 3; tri++) {
                            for (var tri_vtx: u32 = 0; tri_vtx < 3; tri_vtx++) { /* unroll */
                              if (vtx == triangle_indices[tri * 3 + tri_vtx]) {
                                vertex_normals[vtx] += facet_normals[tri];
                              }
                            }
                          }
                          vertex_normals[vtx] = normalize(vertex_normals[vtx]);
                        }
                    }
                  `,
      });