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;  
  int wallFlagged = 0;   
  int nearFlagged = 0;   
  Route route = kTwoPass;
  bool fellBack = false;  
};
 
struct BuildReport {
  long buried = 0;             
  long reachExceeded = 0;      
  long empty = 0;              
  long overflow = 0;           
  long incomplete = 0;         
  int overBufferRebuilds = 0;  
};
 
template <class Real, int MAXP = kMaxPlanes, int MAXT = kMaxTriangles, bool Weighted = false,
          class Sdf = NoSdf>
struct MovingTessellation {
  using Mem = peclet::core::MemSpace;
  using Exec = peclet::core::ExecSpace;
  using PlanePolicy = std::conditional_t<Weighted, Power, Voronoi>;  // radical vs bisector planes
  static constexpr bool kHasSdf = !std::is_same_v<Sdf, NoSdf>;
  // 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 = kSearchWindow, 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 profile = false;  
  int lastOverBufferRebuilds = 0;  
  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;  
 
  Kokkos::View<int*, Mem> lastStatus;  
  BuildReport report() const {
    BuildReport r;
    if (lastStatus.extent(0) == 0)
      return r;
    auto h = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, lastStatus);
    for (int i = 0; i < nProc; ++i) {
      const int s = h(i);
      r.buried += (s & kBuried) ? 1 : 0;
      r.reachExceeded += (s & kReachExceeded) ? 1 : 0;
      r.empty += (s & kEmpty) ? 1 : 0;
      r.overflow += (s & kOverflow) ? 1 : 0;
      r.incomplete += (s & kIncomplete) ? 1 : 0;
    }
    r.overBufferRebuilds = lastOverBufferRebuilds;
    return r;
  }
 
  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; }
 
  // ---- near-miss (face-gain) completeness ----
  // The seed-local certificate cannot see a face GAINED from a seed that is not in the stored
  // topology (its bisector drifts into the cell without either seed moving past the Verlet
  // skin): measured ~0.1 % of neighbour relations missed per step, ~1e-3 volume error over long
  // runs. Every (re)build records the candidates whose plane missed the cell by less than
  // `nearMargin` (a multiple of the skin); the certificate re-tests those planes on the current
  // positions and flags both cells when one now cuts. Set useNearMiss=false for the old (cheaper,
  // incomplete) certificate.
  bool useNearMiss = true;
  int kNearCap = 32;                
  Real nearMarginFrac = Real(0.5);  
  bool nearFresh = false;  
  Kokkos::View<int*, Mem> near, nearCnt;
  Kokkos::View<int, Mem> nearCounter;
 
  // ---- SDF solid (rung A0) ----
  Sdf sdf;                
  WallStore<Real> wall;   
  bool wallExact = true;  
  Real wallSkin = 0;      
  Kokkos::View<int, Mem> wallCounter;
  // 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_ = kSearchWindow,
             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;
    if (useNearMiss) {
      near = Kokkos::View<int*, Mem>(view_alloc(std::string("mt.near"), WithoutInitializing),
                                     (size_t)n * kNearCap);
      nearCnt = Kokkos::View<int*, Mem>("mt.nearCnt", n);
      nearCounter = Kokkos::View<int, Mem>("mt.nearCounter");
    }
    if constexpr (kHasSdf) {
      wall.alloc(n);
      wallCounter = Kokkos::View<int, Mem>("mt.wallCounter");
      surgical = false;  // the surgical re-clip knows nothing of walls; the gather does
    }
  }
 
  KOKKOS_INLINE_FUNCTION static void dispFromRef(const Kokkos::View<Real*, Mem>& P,
                                                 const Kokkos::View<Real*, Mem>& XR, int i, Real Lx,
                                                 Real Ly, Real Lz, Real& dx, Real& dy, Real& dz) {
    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);
    const Real Lxh = Real(0.5) * Lx, Lyh = Real(0.5) * Ly, Lzh = Real(0.5) * Lz;
    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);
  }
 
  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
    // A gate-rebuild (eagerAdj=false) stays in the high-churn regime and skips the near-miss
    // emission like it skips the adjacency; the lists are rebuilt once churn drops (step()).
    const bool emitNear = useNearMiss && eagerAdj;
    auto res = buildTessellation<Real, Weighted, Sdf>(
        pos, weight, N, L, sw, densityCount, gd, sdf, /*withForceGeom=*/false, nBuild, store.np,
        store.nt, store.pnbr, store.tri, /*outCand=*/{}, /*outCandCnt=*/{}, /*candCap=*/0,
        /*wlc=*/nullptr, wall, /*withAreaGrad=*/false, emitNear ? near : Kokkos::View<int*, Mem>{},
        emitNear ? nearCnt : Kokkos::View<int*, Mem>{}, emitNear ? kNearCap : 0,
        nearMarginFrac * skin, /*withWallFD=*/false, /*withMoments=*/false, profile);
    nearFresh = emitNear;
    lastStatus = res.status;
    lastOverBufferRebuilds = res.overBufferRebuilds;
    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>)
    // near-miss (face-gain) check captures
    auto NR = near;
    auto NC = nearCnt;
    auto NCtr = nearCounter;
    const int nearCapL = kNearCap;  // a local: a member read inside the lambda would capture `this`
    const bool nearOn = useNearMiss && useSkin && nearFresh;  // once per step (Pass-1 certify)
    if (nearOn)
      Kokkos::deep_copy(nearCounter, 0);
    // SDF boundary watch captures (force-captured outside the constexpr-if, nvcc rule)
    WallStore<Real> Wl = wall;
    Sdf Sd = sdf;
    auto WC = wallCounter;
    const bool wallExactL = wallExact;
    const Real wallSkin2 = wallSkin * wallSkin;
    if constexpr (kHasSdf)
      Kokkos::deep_copy(wallCounter, 0);
    Kokkos::parallel_for(
        "mt.certify", Kokkos::RangePolicy<Exec>(0, nProc), KOKKOS_LAMBDA(int i) {
          CertCell c;
          st.load(i, c, Lx, Ly, Lz);
          Real dx, dy, dz;
          dispFromRef(P, XR, i, Lx, Ly, Lz, dx, dy, dz);
          const bool wallOk = Wl.load(i, c, dx, dy, dz);  // no-op (true) without a wall store
          // force-capture everything the constexpr boundary watch reads OUTSIDE it (nvcc rule)
          const Sdf SdC = Sd;
          int* wcPtr = WC.data();
          const bool wallExactC = wallExactL;
          const Real wallSkin2C = wallSkin2;
          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) {
            mover = (dx * dx + dy * dy + dz * dz) > half2;
            if (mover)
              Mv(i) = 1;
          }
          // ---- SDF boundary watch (rung A0): what the seed-seed certificate cannot see ----
          bool wflag = false;
          if constexpr (kHasSdf) {
            const Real seed[3] = {P(3 * i), P(3 * i + 1), P(3 * i + 2)};
            if (!wallOk) {
              wflag = true;  // crossed one of its own wall planes
            } else if (Wl.cnt(i) > 0) {
              if (wallExactC || (dx * dx + dy * dy + dz * dz) > wallSkin2C) {
                wflag = true;  // a moved wall cell: only a fresh clip reproduces the cold build
              } else {         // skin mode: keep the planes unless the seed or a vertex entered the
                        // solid (a sagitta-placed plane sits beyond the surface, so the seed can
                        // enter the solid without crossing its own plane)
                if (SdC.eval(seed[0], seed[1], seed[2]) <= Real(0))
                  wflag = true;
                for (int t = 0; t < c.nt && !wflag; ++t)
                  if (c.alive[t] && SdC.eval(seed[0] + c.vx[t], seed[1] + c.vy[t],
                                             seed[2] + c.vz[t]) < Real(-1e-8))
                    wflag = true;
              }
            } else if (c.nt == 0) {
              wflag = SdC.eval(seed[0], seed[1], seed[2]) > Real(0);  // seed back in the fluid
            } else {
              wflag = sdfWouldClip(c, seed, SdC);  // the cold build would cut this cell
            }
            if (wflag)
              Kokkos::atomic_inc(wcPtr);
          }
          // near-miss check: does any recorded non-neighbour plane now cut the cell?
          bool nflag = false;
          if (nearOn) {
            const int nc = NC(i);
            if (nc > nearCapL) {
              nflag = true;  // list overflowed at build: re-gather conservatively
            } else {
              const Real sx = P(3 * i), sy = P(3 * i + 1), sz = P(3 * i + 2);
              for (int q = 0; q < nc; ++q) {
                const int j = NR((size_t)i * nearCapL + q);
                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 pv[3] = {rx, ry, rz};
                Real wN = Real(0);
                if constexpr (Weighted)
                  wN = Wt(j);
                const Real off = PlanePolicy::template offsetFromRel<Real>(pv, wSelfI, wN);
                if (c.planeGap(pv, off) > Real(0)) {
                  nflag = true;
                  if (j < nP_)
                    Kokkos::atomic_exchange(&M(j), 1);  // the gain is mutual
                }
              }
            }
            if (nflag)
              Kokkos::atomic_inc(&NCtr());
          }
          if (!ok || mover || wflag || nflag)
            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;
    WallStore<Real> Wl = wall;
    auto XR = xRef;
    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);
          {
            Real dx, dy, dz;
            dispFromRef(P, XR, i, Lx, Ly, Lz, dx, dy, dz);
            Wl.load(i, c, dx, dy, dz);  // no-op without a wall store
          }
          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;
    WallStore<Real> Wl = wall;
    auto XR = xRef;
    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
          {
            Real dx, dy, dz;
            dispFromRef(P, XR, i, Lx, Ly, Lz, dx, dy, dz);
            Wl.load(i, c, dx, dy, dz);  // no-op without a wall store
          }
          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, Sdf>(grid, wl, n, store.np, store.nt, store.pnbr, store.tri,
                                            vol, store.poke4, sdf, /*withForceGeom=*/false, wall,
                                            near, nearCnt, useNearMiss ? kNearCap : 0,
                                            nearMarginFrac * skin);
    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;
    if constexpr (kHasSdf)
      Kokkos::deep_copy(s.wallFlagged, wallCounter);
    if (useNearMiss)
      Kokkos::deep_copy(s.nearFlagged, nearCounter);
 
    // 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 (useNearMiss && !nearFresh) {
      // First low-churn step after gate-rebuilds: the certificate above ran WITHOUT the near-miss
      // lists (incomplete for face gains), so do not repair on it — pay one emitting rebuild
      // (exact) and repair from the next step on with fresh lists.
      rebuild(pos, /*eagerAdj=*/true);
      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 && !kHasSdf) {
      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