numerics
Loading...
Searching...
No Matches
cell_list.hpp File Reference

Cache-coherent 2D cell list for O(1) amortized neighbour queries. More...

#include <vector>
#include <cmath>
#include <algorithm>
#include <cassert>
#include <utility>

Go to the source code of this file.

Classes

struct  num::IntRange
 Lightweight read-only range over a contiguous int array (C++17-safe span). More...
 
class  num::CellList2D< Scalar >
 

Namespaces

namespace  num
 

Detailed Description

Cache-coherent 2D cell list for O(1) amortized neighbour queries.

Algorithm (CS theory: counting sort + prefix sums)

The standard spatial hash maps cell -> linked list of particles. Each linked-list traversal chases pointers through scattered heap memory, causing a cache miss per particle in the inner neighbour loop.

The cell list replaces this with a counting sort:

build() – O(n + C) where C = total grid cells

  1. Assign each particle i to cell id_i = cy_i*nx + cx_i O(n)
  2. count[id]++ for each particle O(n)
  3. Prefix sum: start[c] = Sigma count[0..c-1] O(C)
  4. Scatter: sorted[start[id_i] + offset++] = i O(n) -> sorted[] contains particle indices grouped by cell. start[c]..start[c+1] is the contiguous range for cell c.

query(px, py, f) – O(k) sequential reads For each of the 9 cells in the 3x3 neighbourhood: iterate sorted[start[c]..start[c+1]) <- sequential memory No pointer chasing. Particles in the same cell sit next to each other in sorted[], so the processor prefetches them correctly.

iterate_pairs(f) – O(n*k/2) Newton's 3rd law half-shell Visits each unique unordered pair {i,j} exactly once. Based on the "forward half-shell" technique from molecular dynamics (Plimpton 1995, LAMMPS):

for each cell (cx, cy): intra-cell: pairs (a, b) where a < b in sorted[] index inter-cell: cross-product with 4 "forward" neighbours: (cx+1, cy-1), (cx+1, cy), (cx+1, cy+1), (cx, cy+1)

These 4 offsets + self cover all 9 neighbour directions uniquely:

  • dx > 0 -> handled when processing the left cell ✓
  • dx = 0, dy > 0 -> handled by (0,+1) offset ✓
  • dx < 0 or dy < 0 -> handled from the other side ✓ Every pair appears in exactly one (cell, forward-offset) combination.
Caller responsibilities
  • Particles must lie within the domain passed to the constructor.
  • get_pos(i) callable must return {x, y} for particle i.
  • Distance check (r < cutoff) must still be done in the callback.
Template Parameters
Scalarfloat or double

Definition in file cell_list.hpp.