     }
   }
 
+  // Build 3 top levels (=apex), which are tiny.
   __global__ __launch_bounds__(128) void reduceApex(uint4* __restrict__ apex_d,
                                                     uint32_t* sum_d,
                                                     const uint32_t* in_d,

    return n;
  }

  template<class FieldType>
  __global__ __launch_bounds__(32) void reduceBaseIndexed(uint8_t* __restrict__           index_cases_d,
                                                          uint4* __restrict__             out_vertex_level0_d,

        float2 next = fetch2(ptr, field_row_stride, cell.y, grid_max_index.y, threshold);
        ptr += cell.z + 1 < grid_max_index.z ? field_slice_stride : 0;

        // Merge bits from previous and current slice.
        float2 t0 = make_float2(__fmaf_rn(float(1 << 4), next.x, prev.x),
                                __fmaf_rn(float(1 << 4), next.y, prev.y));

        prev = next;

        // Fetch case from x+1 using warp-shuffle. Last thread has garbage data,
        // but that result is masked out.
        float2 tt = make_float2(__shfl_down_sync(0xffffffff, t0.x, 1),
                                __shfl_down_sync(0xffffffff, t0.y, 1));

        if (cell.x <= grid_max_index.x && thread < 31) {
          // Merge in results from the warp shuffle. 
          uint32_t g0 = static_cast<uint32_t>(__fmaf_rn(2.f, tt.x, t0.x));
          uint32_t g1 = static_cast<uint32_t>(__fmaf_rn(2.f, tt.y, t0.y));
          uint32_t s0 = __byte_perm(g0, g1, (4 << 12) | (2 << 8) | (1 << 4) | (0 << 0));
          g0 = g0 | (s0 >> 6);
          g1 = g1 | (g1 >> 6);

          // Extract cases for the 5-element y-column of cases that each thread has
          uint32_t case_y0 = __byte_perm(g0, 0, (4 << 12) | (4 << 8) | (4 << 4) | (0 << 0));
          uint32_t case_y1 = __byte_perm(g0, 0, (4 << 12) | (4 << 8) | (4 << 4) | (1 << 0));
          uint32_t case_y2 = __byte_perm(g0, 0, (4 << 12) | (4 << 8) | (4 << 4) | (2 << 0));
          uint32_t case_y3 = __byte_perm(g1, 0, (4 << 12) | (4 << 8) | (4 << 4) | (0 << 0));
          uint32_t case_y4 = __byte_perm(g1, 0, (4 << 12) | (4 << 8) | (4 << 4) | (1 << 0));

          // Look up index count. Set count for cells outside of domain to zero. Careful masking.
          uint32_t ic_y0 = (cell.x < grid_max_index.x) && (cell.z < grid_max_index.z) && ((cell.y + 0u) < grid_max_index.y) ? index_count[case_y0] : 0u;
          uint32_t ic_y1 = (cell.x < grid_max_index.x) && (cell.z < grid_max_index.z) && ((cell.y + 1u) < grid_max_index.y) ? index_count[case_y1] : 0u;
          uint32_t ic_y2 = (cell.x < grid_max_index.x) && (cell.z < grid_max_index.z) && ((cell.y + 2u) < grid_max_index.y) ? index_count[case_y2] : 0u;
          uint32_t ic_y3 = (cell.x < grid_max_index.x) && (cell.z < grid_max_index.z) && ((cell.y + 3u) < grid_max_index.y) ? index_count[case_y3] : 0u;
          uint32_t ic_y4 = (cell.x < grid_max_index.x) && (cell.z < grid_max_index.z) && ((cell.y + 4u) < grid_max_index.y) ? index_count[case_y4] : 0u;

          // Look up vertex count. Set count for cells outside of domain to zero.
          uint32_t vc_y0 = axisCountFromCase(case_y0);
          uint32_t vc_y1 = cell.y + 1 <= grid_max_index.y ? axisCountFromCase(case_y1) : 0;
          uint32_t vc_y2 = cell.y + 2 <= grid_max_index.y ? axisCountFromCase(case_y2) : 0;
