repair.hpp

Go to the documentation of this file.

 
#ifndef PECLET_VORO_REPAIR_HPP
#define PECLET_VORO_REPAIR_HPP
 
#include <cmath>
#include <Kokkos_Core.hpp>
#include <string>
 
#include "peclet/core/common/view.hpp"
#include "peclet/voro/convex_cell.hpp"
#include "peclet/voro/subset_gather.hpp"
#include "peclet/voro/tess_grid.hpp"
#include "peclet/voro/tessellator.hpp"
#include "peclet/voro/topology_store.hpp"
 
namespace peclet::voro {
 
struct RepairStats {
  enum Route { kTwoPass = 0, kDilated = 1, kRebuildGate = 2 };
  int pass1 = 0;  
  int pass1Raw = 0;  
  int pass2 = 0;     
  int extra = 0;     
  int surgical = 0;  
  int verifyPasses = 0;  
  Route route = kTwoPass;
  bool fellBack = false;  
};
 
template <class Real, int MAXP = 64, int MAXT = 112, bool Weighted = false>
struct MovingTessellation {
  using Mem = peclet::core::MemSpace;
  using Exec = peclet::core::ExecSpace;
  using PlanePolicy = std::conditional_t<Weighted, Power, Voronoi>;  // radical vs bisector planes
  // The certificate/load/re-eval path reads the resident `poke4` plane set by pointer and never
  // touches the edge adjacency, so it uses a LEAN cell with NO `adj` member (TrackAdj=false) —
  // ~MAXT*3*4 bytes smaller per thread, which lifts occupancy on the full-N certify (the GPU
  // repair's hot kernel). Only the cell-BUILDING paths (the gather's CellBuilder, rebuildAdjAll,
  // surgicalRepair) need `adj` to derive poke4, so they use the TrackAdj=true BuildCell.
  using CertCell = ConvexCell<Real, MAXP, MAXT, false>;
  using BuildCell = ConvexCell<Real, MAXP, MAXT, true>;
  using Store = TopologyStore<MAXP, MAXT>;
 
  int N = 0;  
  int nProc =
      0;  
  int sw = 4, densityCount = -1;
  Real L[3] = {1, 1, 1};
  Real tol =
      0;  
  Real skin = 0;  
  int verifyCap = 2;  
 
  // ---- Phase 3: adaptive three-way gate (decided after the free certificate, before any gather)
  // ----
  bool useGate = true;  
  bool useDilation =
      false;  
  double churnThresh =
      0.55;  
  double dilateMaxChurn =
      0.15;  
  int clusterNbhd = 18;  
  // ---- Phase 4: single-face surgical repair (no grid gather; gated, default OFF — see note) ----
  // NEGATIVE RESULT (measured): on this engine surgical re-clip is BOTH slower and not robustly
  // exact, so it stays off. (1) Slower: the cold gather clips closest-first with the
  // security-radius early-out (the cell tightens fast, few live triangles per clip), whereas
  // re-clipping the stored+partner candidates unsorted keeps the cell large and costs more O(#tri)
  // horizon scans — the gather wins. (2) Not exact: a surgically-rebuilt cell that GAINS a
  // neighbour absent from its candidate set stays convex, so the certificate can't see the gain
  // (§1c) and the verify loop can't catch it; the small per-step error then COMPOUNDS in the
  // resident store (maxRelV grew to ~0.5 over a sweep). Making it exact needs the maintained
  // 1-/2-ring candidate set (the ConnectivityArena the plan parks). Kept as the implemented
  // experiment behind this flag; the production path is the Phase-3 gated two-pass gather.
  bool surgical =
      false;  
  static constexpr int kExtraCap = 8;  
  // ---- Lawson local convexity certificate (the cheap O(nt) convexity test) ----
  bool localCert =
      true;  
  bool adjFresh =
      false;  
  bool useLocalNow = false;  
 
  Store store;
  Kokkos::View<Real*, Mem> vol;   // N : cell volumes (the published geometry scalar)
  Kokkos::View<Real*, Mem> xRef;  // 3N : positions at the last (re)build (Verlet reference)
  Kokkos::View<Real*, Mem> weight;  // N : per-seed power weights (empty for Voronoi); DOFs alongside
                                    // positions. Set via setWeights before rebuild/step.
 
  void setWeights(const Kokkos::View<Real*, Mem>& w) { weight = w; }
  // scratch
  Kokkos::View<int*, Mem> mask, mask2, mover, rebuilt, wl1, wl2, wlM;
  Kokkos::View<int*, Mem> cellFlag;            // per linear grid-cell flagged-seed count (dilation)
  Kokkos::View<int*, Mem> extraNbr, extraCnt;  // Phase-4 partner-discovered candidate neighbours
  Kokkos::View<int, Mem> counter;
  int gridNcell = 0;
 
  void alloc(int n, const Real Lbox[3], Real tol_, Real skin_, int sw_ = 4, int densityCount_ = -1,
             int nProc_ = -1) {
    N = n;
    nProc = (nProc_ < 0) ? n : nProc_;
    sw = sw_;
    densityCount = densityCount_ < 0 ? n : densityCount_;
    tol = tol_;
    skin = skin_;
    for (int k = 0; k < 3; ++k)
      L[k] = Lbox[k];
    store.alloc(n);
    store.allocPoke4();  // resident local-cert plane set (computed by the gather, read each step by
                         // the cert)
    using Kokkos::view_alloc;
    using Kokkos::WithoutInitializing;
    vol = Kokkos::View<Real*, Mem>("mt.vol", n);
    xRef = Kokkos::View<Real*, Mem>(view_alloc(std::string("mt.xRef"), WithoutInitializing),
                                    (size_t)n * 3);
    mask = Kokkos::View<int*, Mem>("mt.mask", n);
    mask2 = Kokkos::View<int*, Mem>("mt.mask2", n);
    mover = Kokkos::View<int*, Mem>("mt.mover", n);
    rebuilt = Kokkos::View<int*, Mem>("mt.rebuilt", n);
    wl1 = Kokkos::View<int*, Mem>(view_alloc(std::string("mt.wl1"), WithoutInitializing), n);
    wl2 = Kokkos::View<int*, Mem>(view_alloc(std::string("mt.wl2"), WithoutInitializing), n);
    wlM = Kokkos::View<int*, Mem>(view_alloc(std::string("mt.wlM"), WithoutInitializing), n);
    extraNbr = Kokkos::View<int*, Mem>(view_alloc(std::string("mt.extraNbr"), WithoutInitializing),
                                       (size_t)n * kExtraCap);
    extraCnt = Kokkos::View<int*, Mem>("mt.extraCnt", n);
    counter = Kokkos::View<int, Mem>("mt.counter");
    // Per-backend gate threshold: rebuild is cheap + clip-bound on GPU (repair loses past ~50%
    // flagged), expensive on the CPU paths (repair pays out to ~85%). Hand-set from the
    // cross-device sweep; the plan's Phase-3 tuning would fit these from logged (signals vs
    // verify-clean) per backend.
    constexpr bool host = Kokkos::SpaceAccessibility<Kokkos::HostSpace, Mem>::accessible;
    churnThresh = host ? 0.70 : 0.50;
  }
 
  void rebuild(const Kokkos::View<Real*, Mem>& pos, bool eagerAdj = true) {
    Kokkos::View<long*, Mem> gd;
    const int nBuild =
        (nProc == N) ? -1 : nProc;  // build only the owned cells; ghosts are candidates
    auto res = buildTessellation<Real, Weighted>(pos, weight, N, L, sw, densityCount, gd, NoSdf{},
                                                 /*withForceGeom=*/false, nBuild, store.np, store.nt,
                                                 store.pnbr, store.tri);
    Kokkos::deep_copy(vol, res.view.cellVolume);
    Kokkos::deep_copy(xRef, pos);
    adjFresh = false;  // the cold build does not emit adjacency
    if (eagerAdj)
      maybeBuildAdj(pos);  // the whole topology changed -> repopulate store.poke4 (unless gating)
  }
 
  int compact(const Kokkos::View<int*, Mem>& m, const Kokkos::View<int*, Mem>& wl) {
    Kokkos::deep_copy(counter, 0);
    auto M = m;
    auto W = wl;
    auto c = counter;
    Kokkos::parallel_for(
        "mt.compact", Kokkos::RangePolicy<Exec>(0, nProc), KOKKOS_LAMBDA(int i) {
          if (M(i))
            W(Kokkos::atomic_fetch_add(&c(), 1)) = i;
        });
    int h = 0;
    Kokkos::deep_copy(h, counter);
    return h;
  }
 
  template <bool Local>
  void certify(const Kokkos::View<Real*, Mem>& pos, const Kokkos::View<int*, Mem>& outMask,
               const Kokkos::View<int*, Mem>& outMover, bool useSkin) {
    Kokkos::deep_copy(outMask, 0);
    if (useSkin)
      Kokkos::deep_copy(outMover, 0);
    // Phase-4: in the Pass-1 certify, also record partner-discovered new neighbours for the
    // surgical path. The gaining cells of a flip are the violated-plane partners; each pair gains
    // the other, so scatter every other partner of a flagging cell onto each partner's
    // extra-candidate list (no dedup — duplicates are harmless no-op clips; overflow past kExtraCap
    // just truncates and the cell falls back to a full gather at verify). Only in Pass 1 (useSkin),
    // and only when surgical is enabled.
    const bool doExtra = surgical && useSkin;
    if (doExtra)
      Kokkos::deep_copy(extraCnt, 0);
    auto EN = extraNbr;
    auto EC = extraCnt;
    constexpr int kCap = kExtraCap;
    auto P = pos;
    auto XR = xRef;
    auto Vv = vol;
    auto M = outMask;
    auto Mv = outMover;
    Store st = store;
    auto Pk4 = store.poke4;
    const Real Lx = L[0], Ly = L[1], Lz = L[2];
    const Real Lxh = Real(0.5) * Lx, Lyh = Real(0.5) * Ly, Lzh = Real(0.5) * Lz;
    const Real tolL = tol, half2 = Real(0.25) * skin * skin;
    const bool skinOn = useSkin;
    const int nP_ = nProc;
    auto Wt = weight;  // power weights (empty for Voronoi; captured for reevalGeometry<PlanePolicy>)
    Kokkos::parallel_for(
        "mt.certify", Kokkos::RangePolicy<Exec>(0, nProc), KOKKOS_LAMBDA(int i) {
          CertCell c;
          st.load(i, c, Lx, Ly, Lz);
          const Real* wPtr = Wt.data();  // force-capture Wt OUTSIDE the constexpr-if (nvcc rule)
          Real wSelfI = Real(0);
          if constexpr (Weighted) wSelfI = Wt(i);
          c.template reevalGeometry<PlanePolicy>(P(3 * i), P(3 * i + 1), P(3 * i + 2), P.data(), Lx,
                                                 wSelfI, wPtr);
          Vv(i) = c.volumePerVertex();
          int partners[32];
          int nP = 0;
          bool ok;
          // Force-capture every lambda capture used in the branches OUTSIDE the constexpr-if (nvcc:
          // an extended __host__ __device__ lambda cannot FIRST-capture a variable inside a
          // constexpr-if).
          const unsigned char* pk4 = &Pk4((size_t)i * MAXT * 4);  // poke4 is always allocated
          const Real tl = tolL;
          if constexpr (Local)
            ok = c.isLocallyConvexPartners(pk4, tl, partners, 32, nP);
          else {
            (void)pk4;
            ok = c.isSelfConsistentPartners(tl, partners, 32, nP);
          }
          bool mover = false;
          if (skinOn) {
            Real dx = P(3 * i) - XR(3 * i), dy = P(3 * i + 1) - XR(3 * i + 1),
                 dz = P(3 * i + 2) - XR(3 * i + 2);
            dx = dx > Lxh ? dx - Lx : (dx < -Lxh ? dx + Lx : dx);
            dy = dy > Lyh ? dy - Ly : (dy < -Lyh ? dy + Ly : dy);
            dz = dz > Lzh ? dz - Lz : (dz < -Lzh ? dz + Lz : dz);
            mover = (dx * dx + dy * dy + dz * dz) > half2;
            if (mover)
              Mv(i) = 1;
          }
          if (!ok || mover)
            M(i) = 1;
          // only mark partners we MAINTAIN (owned); a ghost partner is the owning rank's
          // responsibility.
          for (int q = 0; q < nP; ++q)
            if (partners[q] < nP_)
              Kokkos::atomic_exchange(&M(partners[q]), 1);
          if (doExtra) {  // each pair of gaining cells gains the other -> seed the surgical
                          // candidate lists
            for (int q = 0; q < nP; ++q) {
              const int cell = partners[q];
              if (cell < 0 || cell >= nP_)
                continue;
              for (int r = 0; r < nP; ++r) {
                if (r == q)
                  continue;
                const int nb = partners[r];
                if (nb < 0)
                  continue;
                const int slot = Kokkos::atomic_fetch_add(&EC(cell), 1);
                if (slot < kCap)
                  EN((size_t)cell * kCap + slot) = nb;
              }
            }
          }
        });
  }
 
  template <bool Local>
  void certifyList(const Kokkos::View<Real*, Mem>& pos, const Kokkos::View<int*, Mem>& list, int n,
                   const Kokkos::View<int*, Mem>& outMask) {
    Kokkos::deep_copy(outMask, 0);
    if (n <= 0)
      return;
    auto P = pos;
    auto Vv = vol;
    auto M = outMask;
    Store st = store;
    auto Pk4 = store.poke4;
    const Real Lx = L[0], Ly = L[1], Lz = L[2];
    const Real tolL = tol;
    const int nP_ = nProc;
    auto Wt = weight;
    Kokkos::parallel_for(
        "mt.certifyList", Kokkos::RangePolicy<Exec>(0, n), KOKKOS_LAMBDA(int s) {
          const int i = list(s);
          CertCell c;
          st.load(i, c, Lx, Ly, Lz);
          const Real* wPtr = Wt.data();  // force-capture Wt OUTSIDE the constexpr-if (nvcc rule)
          Real wSelfI = Real(0);
          if constexpr (Weighted) wSelfI = Wt(i);
          c.template reevalGeometry<PlanePolicy>(P(3 * i), P(3 * i + 1), P(3 * i + 2), P.data(), Lx,
                                                 wSelfI, wPtr);
          Vv(i) = c.volumePerVertex();
          int partners[32];
          int nP = 0;
          bool ok;
          const unsigned char* pk4 =
              &Pk4((size_t)i * MAXT * 4);  // force-capture poke4 + tol OUTSIDE the
          const Real tl = tolL;            // constexpr-if (nvcc extended-lambda rule)
          if constexpr (Local)
            ok = c.isLocallyConvexPartners(pk4, tl, partners, 32, nP);
          else {
            (void)pk4;
            ok = c.isSelfConsistentPartners(tl, partners, 32, nP);
          }
          if (!ok)
            M(i) = 1;
          for (int q = 0; q < nP; ++q)
            if (partners[q] < nP_)
              Kokkos::atomic_exchange(&M(partners[q]), 1);
        });
  }
 
  void rebuildAdjAll(const Kokkos::View<Real*, Mem>& pos) {
    if (!localCert)
      return;
    Store st = store;
    auto P = pos;
    auto Pk4 = store.poke4;
    const Real Lx = L[0], Ly = L[1], Lz = L[2];
    auto Wt = weight;
    Kokkos::parallel_for(
        "mt.rebuildAdjAll", Kokkos::RangePolicy<Exec>(0, nProc), KOKKOS_LAMBDA(int i) {
          BuildCell c;
          st.load(i, c, Lx, Ly,
                  Lz);  // load brings topology; adj is rebuilt to derive the cert planes
          const Real* wPtr = Wt.data();  // force-capture Wt OUTSIDE the constexpr-if (nvcc rule)
          Real wSelfI = Real(0);
          if constexpr (Weighted) wSelfI = Wt(i);
          c.template reevalGeometry<PlanePolicy>(P(3 * i), P(3 * i + 1), P(3 * i + 2), P.data(), Lx,
                                                 wSelfI, wPtr);  // computePoke4 needs vertices
          c.rebuildAdjacency();
          c.computePoke4(&Pk4((size_t)i * MAXT *
                              4));  // the only thing this pass produces (topology unchanged)
        });
  }
  void maybeBuildAdj(const Kokkos::View<Real*, Mem>& pos) {
    if (localCert && !adjFresh) {
      rebuildAdjAll(pos);
      adjFresh = true;
      useLocalNow = true;
    }
  }
 
  void gatherSet(const TessGrid<Real>& grid, const Kokkos::View<int*, Mem>& wl, int n,
                 const Kokkos::View<Real*, Mem>& pos) {
    if (n <= 0)
      return;
    // TrackAdj gather: clip stitches the edge adjacency incrementally; the builder derives the cert
    // plane set from it (computePoke4, no findSharing) into store.poke4, so the local cert stays
    // valid for the gathered cells with no separate maintenance pass.
    subsetGather<Real, Weighted, true>(grid, wl, n, store.np, store.nt, store.pnbr, store.tri, vol,
                                       store.poke4, NoSdf{}, /*withForceGeom=*/false);
    auto W = wl;
    auto rb = rebuilt;
    auto P = pos;
    auto XR = xRef;
    Kokkos::parallel_for(
        "mt.markRebuilt", Kokkos::RangePolicy<Exec>(0, n), KOKKOS_LAMBDA(int s) {
          const int i = W(s);
          rb(i) = 1;
          XR(3 * i) = P(3 * i);
          XR(3 * i + 1) = P(3 * i + 1);
          XR(3 * i + 2) = P(3 * i + 2);
        });
  }
 
  void collectNewNbrs(const Kokkos::View<int*, Mem>& wl, int n,
                      const Kokkos::View<int*, Mem>& outMask) {
    Kokkos::deep_copy(outMask, 0);
    if (n <= 0)
      return;
    Store st = store;
    auto rb = rebuilt;
    auto M = outMask;
    const int nP_ = nProc;
    Kokkos::parallel_for(
        "mt.collectNbrs", Kokkos::RangePolicy<Exec>(0, n), KOKKOS_LAMBDA(int s) {
          const int i = wl(s);
          const int np = st.np(i), nt = st.nt(i);
          for (int k = 6; k < np; ++k) {
            const int j = st.pnbr((size_t)i * MAXP + k);
            if (j < 0 || j >= nP_ || rb(j))
              continue;   // only expand neighbours we MAINTAIN (owned)
            int cnt = 0;  // face = ≥3 incident live triangles
            for (int t = 0; t < nt; ++t) {
              const unsigned w = st.tri((size_t)i * MAXT + t);
              if (!((w >> 24) & 1u))
                continue;
              const int a = (int)(w & 0xffu), b = (int)((w >> 8) & 0xffu),
                        cc = (int)((w >> 16) & 0xffu);
              if (a == k || b == k || cc == k)
                ++cnt;
            }
            if (cnt >= 3)
              Kokkos::atomic_exchange(&M(j), 1);
          }
        });
  }
 
  void surgicalRepair(const Kokkos::View<int*, Mem>& wl, int n,
                      const Kokkos::View<Real*, Mem>& pos) {
    if (n <= 0)
      return;
    auto W = wl;
    auto P = pos;
    auto rb = rebuilt;
    auto XR = xRef;
    auto Vv = vol;
    auto EN = extraNbr;
    auto EC = extraCnt;
    constexpr int kCap = kExtraCap;
    Store st = store;
    auto Pk4 = store.poke4;
    const Real Lx = L[0], Ly = L[1], Lz = L[2];
    const Real Lxh = Real(0.5) * Lx, Lyh = Real(0.5) * Ly, Lzh = Real(0.5) * Lz;
    Kokkos::parallel_for(
        "mt.surgical", Kokkos::RangePolicy<Exec>(0, n), KOKKOS_LAMBDA(int s) {
          const int i = W(s);
          BuildCell c;
          c.initBox(Lx, Ly, Lz);
          const Real sx = P(3 * i), sy = P(3 * i + 1), sz = P(3 * i + 2);
          const int npi = st.np(i);
          for (int k = 6; k < npi && !c.overflow; ++k) {  // stored neighbours, current positions
            const int j = st.pnbr((size_t)i * MAXP + k);
            if (j < 0)
              continue;
            Real rx = P(3 * j) - sx, ry = P(3 * j + 1) - sy, rz = P(3 * j + 2) - sz;
            rx = rx > Lxh ? rx - Lx : (rx < -Lxh ? rx + Lx : rx);
            ry = ry > Lyh ? ry - Ly : (ry < -Lyh ? ry + Ly : ry);
            rz = rz > Lzh ? rz - Lz : (rz < -Lzh ? rz + Lz : rz);
            const Real nrm[3] = {rx, ry, rz};
            c.clip(nrm, Real(0.5) * (rx * rx + ry * ry + rz * rz), j);
          }
          const int ne = EC(i) < kCap ? EC(i) : kCap;
          for (int e = 0; e < ne && !c.overflow; ++e) {  // partner-discovered new neighbours
            const int j = EN((size_t)i * kCap + e);
            if (j < 0 || j == i)
              continue;
            Real rx = P(3 * j) - sx, ry = P(3 * j + 1) - sy, rz = P(3 * j + 2) - sz;
            rx = rx > Lxh ? rx - Lx : (rx < -Lxh ? rx + Lx : rx);
            ry = ry > Lyh ? ry - Ly : (ry < -Lyh ? ry + Ly : ry);
            rz = rz > Lzh ? rz - Lz : (rz < -Lzh ? rz + Lz : rz);
            const Real nrm[3] = {rx, ry, rz};
            c.clip(nrm, Real(0.5) * (rx * rx + ry * ry + rz * rz), j);
          }
          if (!c.overflow) {
            st.save(i, c);
            Vv(i) = c.volumePerVertex();
            c.computePoke4(
                &Pk4((size_t)i * MAXT * 4));  // refresh the cert plane set for this re-clipped cell
          }
          rb(i) = 1;
          XR(3 * i) = sx;
          XR(3 * i + 1) = sy;
          XR(3 * i + 2) = sz;
        });
  }
 
  bool dilate(const TessGrid<Real>& grid, const Kokkos::View<Real*, Mem>& pos) {
    const int dimx = grid.dimx, dimy = grid.dimy, dimz = grid.dimz;
    const int ncell = dimx * dimy * dimz;
    if ((int)cellFlag.extent(0) < ncell)
      cellFlag = Kokkos::View<int*, Mem>(
          Kokkos::view_alloc(std::string("mt.cellFlag"), Kokkos::WithoutInitializing), ncell);
    Kokkos::deep_copy(cellFlag, 0);
    const Real icx = grid.icx, icy = grid.icy, icz = grid.icz;
    auto P = pos;
    auto M = mask;
    auto CF = cellFlag;
    Kokkos::parallel_for(
        "mt.dilate.count", Kokkos::RangePolicy<Exec>(0, nProc), KOKKOS_LAMBDA(int i) {
          if (!M(i))
            return;
          const int gx = ((int)Kokkos::floor(P(3 * i) * icx) % dimx + dimx) % dimx;
          const int gy = ((int)Kokkos::floor(P(3 * i + 1) * icy) % dimy + dimy) % dimy;
          const int gz = ((int)Kokkos::floor(P(3 * i + 2) * icz) % dimz + dimz) % dimz;
          Kokkos::atomic_inc(&CF(gx + gy * dimx + gz * dimx * dimy));
        });
    const int thr = clusterNbhd;
    long added = 0;
    Kokkos::parallel_reduce(
        "mt.dilate.mark", Kokkos::RangePolicy<Exec>(0, nProc),
        KOKKOS_LAMBDA(int i, long& a) {
          if (M(i))
            return;
          const int gx = ((int)Kokkos::floor(P(3 * i) * icx) % dimx + dimx) % dimx;
          const int gy = ((int)Kokkos::floor(P(3 * i + 1) * icy) % dimy + dimy) % dimy;
          const int gz = ((int)Kokkos::floor(P(3 * i + 2) * icz) % dimz + dimz) % dimz;
          int sum = 0;
          for (int dz = -1; dz <= 1; ++dz)
            for (int dy = -1; dy <= 1; ++dy)
              for (int dx = -1; dx <= 1; ++dx) {
                const int nx = (gx + dx + dimx) % dimx, ny = (gy + dy + dimy) % dimy,
                          nz = (gz + dz + dimz) % dimz;
                sum += CF(nx + ny * dimx + nz * dimx * dimy);
              }
          if (sum > thr) {
            M(i) = 1;
            ++a;
          }
        },
        added);
    return added > 0;
  }
 
  void certifyDispatch(const Kokkos::View<Real*, Mem>& pos, const Kokkos::View<int*, Mem>& outMask,
                       const Kokkos::View<int*, Mem>& outMover, bool useSkin) {
    if (useLocalNow)
      certify<true>(pos, outMask, outMover, useSkin);
    else
      certify<false>(pos, outMask, outMover, useSkin);
  }
  void certifyListDispatch(const Kokkos::View<Real*, Mem>& pos, const Kokkos::View<int*, Mem>& list,
                           int n, const Kokkos::View<int*, Mem>& outMask) {
    if (useLocalNow)
      certifyList<true>(pos, list, n, outMask);
    else
      certifyList<false>(pos, list, n, outMask);
  }
 
  RepairStats step(const Kokkos::View<Real*, Mem>& pos) {
    RepairStats s;
    auto grid = buildTessGrid<Real, Weighted>(pos, weight, N, L, sw, densityCount,
                                              Kokkos::View<long*, Mem>());
    Kokkos::deep_copy(rebuilt, 0);
    // Use the cheap local certificate whenever store.poke4 is valid for the current topology. With
    // the 4th-face plane included the local cert matches the brute flag set (it is complete, not a
    // subset), so no per-step displacement gate or brute fallback is needed — large moves and
    // teleports are handled by the same detect→gather→verify path. (adj/face4 are invalidated only
    // by a cold rebuild.)
    useLocalNow = localCert && adjFresh;
 
    // Pass 1: detect (flagged ∪ partners ∪ skin-movers). flagged-commons + their violated-plane
    // partners (the gaining cells) fully cover an isolated flip — no Pass 2 for flips.
    certifyDispatch(pos, mask, mover, /*useSkin=*/true);
    int n1 = compact(mask, wl1);
    s.pass1Raw = n1;
 
    // Small-displacement fast path: the certificate found every cell consistent and no skin-mover
    // (n1 counts flagged ∪ partners ∪ movers). By the §1 completeness argument a stale topology
    // always flags some cell, so n1==0 ⇒ the stored topology is still Voronoi on the new positions
    // — the re-eval already refreshed every volume, so publish with NO gather and NO verify pass.
    // This is the "almost no cells invalidated" regime: the per-step cost collapses to one re-eval
    // + the grid.
    if (n1 == 0) {
      maybeBuildAdj(pos);
      return s;
    }
 
    // Phase-3 adaptive gate (decided here, after the FREE certificate, before any gather):
    //   high global churn  -> a full rebuild is cheaper than repairing this many cells;
    //   dense local cluster -> dilate the Pass-1 set by a grid-cell buffer (one regional gather);
    //   else (sparse)       -> the two-pass repair below.
    if (useGate && (double)n1 / (nProc > 0 ? nProc : 1) > churnThresh) {
      rebuild(pos,
              /*eagerAdj=*/false);  // high-churn regime: stay on brute, don't pay adj maintenance
      s.route = RepairStats::kRebuildGate;
      return s;
    }
    if (useGate && useDilation && n1 > 0 && (double)n1 / nProc < dilateMaxChurn &&
        dilate(grid, pos)) {
      n1 = compact(mask, wl1);
      s.route = RepairStats::kDilated;
    }
    s.pass1 = n1;
    int nM = compact(mover, wlM);
    // Pass 1 gather. Phase 4: FLIP cells (flagged non-movers) are repaired SURGICALLY (re-clip
    // known candidates, no grid gather); far-mover insertions keep the full gather (their new
    // neighbours are in no candidate list). Default (surgical off): the whole Pass-1 set goes
    // through the full gather.
    if (surgical) {
      auto Msk = mask;
      auto Mv = mover;
      auto M2 = mask2;
      Kokkos::parallel_for(
          "mt.surgMask", Kokkos::RangePolicy<Exec>(0, nProc),
          KOKKOS_LAMBDA(int i) { M2(i) = (Msk(i) && !Mv(i)) ? 1 : 0; });
      int nS =
          compact(mask2, wl1);  // wl1's Pass-1 list no longer needed -> reuse for the flip list
      s.surgical = nS;
      surgicalRepair(wl1, nS, pos);   // flips: no gather
      gatherSet(grid, wlM, nM, pos);  // movers: full gather
    } else {
      gatherSet(grid, wl1, n1, pos);
    }
 
    // Pass 2: ONLY the new face-neighbours of the skin-MOVERS (the insertion side §1c — those
    // neighbours never flag and are never partners). Expanding every Pass-1 cell would re-clip ~all
    // neighbours of every flip cluster for nothing; movers are the only cells whose rebuild reveals
    // an unflagged gainer.
    collectNewNbrs(wlM, nM, mask2);
    int n2 = compact(mask2, wl2);
    s.pass2 = n2;
    gatherSet(grid, wl2, n2, pos);
    maybeBuildAdj(
        pos);  // first low-churn step after a rebuild: populate adj so verify + next steps go local
 
    // Verify (scoped): re-certify only the GATHERED cells (every non-gathered cell is provably
    // unchanged — see certifyList). Clean ⇒ publish; else gather the still-flagged ∪ partners
    // (residual coupled flips / tol-degenerate cells) and re-verify, bounded.
    for (int vp = 0; vp < verifyCap + 1; ++vp) {
      const int ng = compact(rebuilt, wl1);  // all cells gathered so far this step
      // Verify with the step's certificate: a far/teleport step already forced useLocalNow=false,
      // so its verify is the COMPLETE brute form (catches far-pokes the gather left on a gathered
      // cell). A small- displacement local step uses the cheap local verify over the
      // freshly-gathered (locally convex) set.
      certifyListDispatch(pos, wl1, ng, mask);
      int nv = compact(mask, wl2);
      s.verifyPasses = vp + 1;
      if (nv == 0)
        return s;  // clean — done
      if (vp == verifyCap)
        break;  // out of budget — fall back
      s.extra += nv;
      gatherSet(grid, wl2, nv, pos);
    }
 
    // Still dirty after the verify budget: cold rebuild (always exact). Rare; the data point that
    // says "for this step, rebuild was cheaper" — the adaptive gate (Phase 3) will pre-empt it up
    // front.
    rebuild(pos);
    s.fellBack = true;
    return s;
  }
};
 
}  // namespace peclet::voro
 
#endif  // PECLET_VORO_REPAIR_HPP