convex_cell.hpp

Go to the documentation of this file.

 
#ifndef PECLET_VORO_CONVEX_CELL_HPP
#define PECLET_VORO_CONVEX_CELL_HPP
 
#include <Kokkos_Core.hpp>
 
#include "peclet/voro/plane_policy.hpp"  // Voronoi/Power plane-from-DOF policies (leaf: Kokkos only)
 
// Loop unroll hint, applied ONLY in the CUDA/HIP device passes (where it lets the compiler
// scalarize small dynamically-indexed local arrays into registers); a no-op on the host so gcc sees
// no unknown pragma.
#if defined(__CUDA_ARCH__) || defined(__HIP_DEVICE_COMPILE__)
#define VOR_UNROLL _Pragma("unroll")
#else
#define VOR_UNROLL
#endif
 
namespace peclet::voro {
 
namespace detail {
// Minimal forward-mode dual number (value + K partials) used by ConvexCell::geomVolumeAreaGrad (the
// geomFull area-Jacobian). Templated on Real so it routes through nvcc/hipcc like the rest of the
// header.
template <class Real, int K>
struct Dual {
  Real v;
  Real d[K];
};
template <class Real, int K>
KOKKOS_INLINE_FUNCTION Dual<Real, K> dnum(Real c) {
  Dual<Real, K> r;
  r.v = c;
  for (int i = 0; i < K; ++i)
    r.d[i] = Real(0);
  return r;
}
template <class Real, int K>
KOKKOS_INLINE_FUNCTION Dual<Real, K> dseed(Real c, int s) {
  Dual<Real, K> r = dnum<Real, K>(c);
  r.d[s] = Real(1);
  return r;
}
template <class Real, int K>
KOKKOS_INLINE_FUNCTION Dual<Real, K> operator+(const Dual<Real, K>& a, const Dual<Real, K>& b) {
  Dual<Real, K> r;
  r.v = a.v + b.v;
  for (int i = 0; i < K; ++i)
    r.d[i] = a.d[i] + b.d[i];
  return r;
}
template <class Real, int K>
KOKKOS_INLINE_FUNCTION Dual<Real, K> operator-(const Dual<Real, K>& a, const Dual<Real, K>& b) {
  Dual<Real, K> r;
  r.v = a.v - b.v;
  for (int i = 0; i < K; ++i)
    r.d[i] = a.d[i] - b.d[i];
  return r;
}
template <class Real, int K>
KOKKOS_INLINE_FUNCTION Dual<Real, K> operator*(const Dual<Real, K>& a, const Dual<Real, K>& b) {
  Dual<Real, K> r;
  r.v = a.v * b.v;
  for (int i = 0; i < K; ++i)
    r.d[i] = a.d[i] * b.v + a.v * b.d[i];
  return r;
}
template <class Real, int K>
KOKKOS_INLINE_FUNCTION Dual<Real, K> operator/(const Dual<Real, K>& a, const Dual<Real, K>& b) {
  Dual<Real, K> r;
  const Real inv = Real(1) / b.v;
  r.v = a.v * inv;
  for (int i = 0; i < K; ++i)
    r.d[i] = (a.d[i] - r.v * b.d[i]) * inv;
  return r;
}
template <class Real, int K>
KOKKOS_INLINE_FUNCTION Dual<Real, K> dsqrt(const Dual<Real, K>& a) {
  Dual<Real, K> r;
  r.v = Kokkos::sqrt(a.v);
  const Real inv = (r.v > Real(0)) ? Real(1) / (Real(2) * r.v) : Real(0);
  for (int i = 0; i < K; ++i)
    r.d[i] = a.d[i] * inv;
  return r;
}
template <class Real, int K>
KOKKOS_INLINE_FUNCTION Dual<Real, K> ddot(const Dual<Real, K> a[3], const Dual<Real, K> b[3]) {
  return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
}
template <class Real, int K>
KOKKOS_INLINE_FUNCTION void dcross(const Dual<Real, K> a[3], const Dual<Real, K> b[3],
                                   Dual<Real, K> o[3]) {
  o[0] = a[1] * b[2] - a[2] * b[1];
  o[1] = a[2] * b[0] - a[0] * b[2];
  o[2] = a[0] * b[1] - a[1] * b[0];
}
template <class Real, int K>
KOKKOS_INLINE_FUNCTION Dual<Real, K> ddet3(const Dual<Real, K> a[3], const Dual<Real, K> b[3],
                                           const Dual<Real, K> c[3]) {
  Dual<Real, K> bc[3];
  dcross(b, c, bc);
  return ddot(a, bc);
}
template <class Real, int K>
KOKKOS_INLINE_FUNCTION void dedgeFoot(const Dual<Real, K> v[3], const Dual<Real, K> c[3],
                                      Dual<Real, K> f[3]) {
  const Dual<Real, K> s = ddot(v, c) / ddot(c, c);
  f[0] = v[0] - s * c[0];
  f[1] = v[1] - s * c[1];
  f[2] = v[2] - s * c[2];
}
}  // namespace detail
 
template <class Real, int MAXP = 64, int MAXT = 96, bool TrackAdj = false>
struct ConvexCell {
  static_assert(MAXP <= 255, "plane index must fit in unsigned char");
  static constexpr bool kTrackAdj =
      TrackAdj;  
  Real n[MAXP]
        [3];      
  Real nn[MAXP];  
  int pnbr[MAXP];                              
  int np;                                      
  unsigned char t0[MAXT], t1[MAXT], t2[MAXT];  
  Real vx[MAXT], vy[MAXT], vz[MAXT];           
  bool alive[MAXT];                            
  int nt;                                      
  bool overflow;  
 
  Kokkos::Array<int, (TrackAdj ? MAXT : 0) * 3> adj;
 
  KOKKOS_INLINE_FUNCTION void initBox(Real L0, Real L1, Real L2) {
    const Real h[3] = {Real(0.5) * L0, Real(0.5) * L1, Real(0.5) * L2};
    for (int ax = 0; ax < 3; ++ax) {
      for (int s = 0; s < 2; ++s) {
        const int k = 2 * ax + s;
        n[k][0] = n[k][1] = n[k][2] = 0;
        n[k][ax] = (s == 0) ? h[ax] : -h[ax];  // foot point = (±h) e_ax
        nn[k] = h[ax] * h[ax];
        pnbr[k] = -1;
      }
    }
    np = 6;
    nt = 0;
    overflow = false;
    // 8 corners: x-plane (0|1), y-plane (2|3), z-plane (4|5)
    for (int sx = 0; sx < 2; ++sx)
      for (int sy = 0; sy < 2; ++sy)
        for (int sz = 0; sz < 2; ++sz) {
          t0[nt] = (unsigned char)(0 + sx);
          t1[nt] = (unsigned char)(2 + sy);
          t2[nt] = (unsigned char)(4 + sz);
          alive[nt] = true;
          computeVertex(nt);
          ++nt;
        }
    if constexpr (TrackAdj)
      rebuildAdjacency();  // one-time bootstrap (~12 box triangles)
  }
 
  KOKKOS_INLINE_FUNCTION void computeVertex(int t) {
    const Real* a = n[t0[t]];
    const Real* b = n[t1[t]];
    const Real* c = n[t2[t]];
    const Real da = nn[t0[t]], db = nn[t1[t]], dc = nn[t2[t]];
    const Real bc[3] = {b[1] * c[2] - b[2] * c[1], b[2] * c[0] - b[0] * c[2],
                        b[0] * c[1] - b[1] * c[0]};
    const Real ca[3] = {c[1] * a[2] - c[2] * a[1], c[2] * a[0] - c[0] * a[2],
                        c[0] * a[1] - c[1] * a[0]};
    const Real ab[3] = {a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2],
                        a[0] * b[1] - a[1] * b[0]};
    const Real det = a[0] * bc[0] + a[1] * bc[1] + a[2] * bc[2];
    const Real inv = det != Real(0) ? Real(1) / det : Real(0);
    vx[t] = (da * bc[0] + db * ca[0] + dc * ab[0]) * inv;
    vy[t] = (da * bc[1] + db * ca[1] + dc * ab[1]) * inv;
    vz[t] = (da * bc[2] + db * ca[2] + dc * ab[2]) * inv;
  }
 
  KOKKOS_INLINE_FUNCTION void initBoxPlanes(Real L0, Real L1, Real L2) {
    const Real h[3] = {Real(0.5) * L0, Real(0.5) * L1, Real(0.5) * L2};
    for (int ax = 0; ax < 3; ++ax)
      for (int s = 0; s < 2; ++s) {
        const int k = 2 * ax + s;
        n[k][0] = n[k][1] = n[k][2] = 0;
        n[k][ax] = (s == 0) ? h[ax] : -h[ax];  // foot point = (±h) e_ax
        nn[k] = h[ax] * h[ax];
        pnbr[k] = -1;
      }
    np = 6;
    overflow = false;
  }
 
  template <class Policy = Voronoi>
  KOKKOS_INLINE_FUNCTION void reevalGeometry(Real sx, Real sy, Real sz, const Real* pos, Real L,
                                             Real wSelf = Real(0), const Real* weight = nullptr) {
    const Real Lh = Real(0.5) * L;
    for (int k = 0; k < np; ++k) {
      if (pnbr[k] < 0)
        continue;  // bounding-box / SDF-wall plane — does not move
      const int g = pnbr[k];
      Real rx = pos[3 * g] - sx, ry = pos[3 * g + 1] - sy, rz = pos[3 * g + 2] - sz;
      rx = rx > Lh ? rx - L : (rx < -Lh ? rx + L : rx);
      ry = ry > Lh ? ry - L : (ry < -Lh ? ry + L : ry);
      rz = rz > Lh ? rz - L : (rz < -Lh ? rz + L : rz);
      if constexpr (Policy::kHasWeightDof) {
        // radical foot point: n = (d/|r|²) r,  d = ½(|r|² + w_self − w_nbr)
        const Real r[3] = {rx, ry, rz};
        const Real rho = rx * rx + ry * ry + rz * rz;
        const Real off = Policy::template offsetFromRel<Real>(r, wSelf, weight[g]);
        const Real a = rho > Real(0) ? off / rho : Real(0);
        n[k][0] = a * rx;
        n[k][1] = a * ry;
        n[k][2] = a * rz;
        nn[k] = n[k][0] * n[k][0] + n[k][1] * n[k][1] + n[k][2] * n[k][2];
      } else {
        // foot point of the Voronoi bisector: n = ½r, offset nn = |n|² = ¼|r|²
        const Real hx = Real(0.5) * rx, hy = Real(0.5) * ry, hz = Real(0.5) * rz;
        n[k][0] = hx;
        n[k][1] = hy;
        n[k][2] = hz;
        nn[k] = hx * hx + hy * hy + hz * hz;
      }
    }
    for (int t = 0; t < nt; ++t)
      if (alive[t])
        computeVertex(t);
  }
 
  KOKKOS_INLINE_FUNCTION bool isSelfConsistent(Real tol) const {
    // Branch-free max-slack accumulation (no per-thread early-out divergence): the whole-cell worst
    // perpendicular poke is compared to `tol` once at the end. The brute scan is the complete
    // certificate — it catches far-pokes by NON-adjacent planes too, which the local (Lawson) form
    // deliberately does not. (An early-out variant would only help when np is large AND fail-fast
    // rate is high; measure first.)
    Real maxd = Real(0);
    for (int t = 0; t < nt; ++t) {
      if (!alive[t])
        continue;
      const Real v[3] = {vx[t], vy[t], vz[t]};
      for (int k = 0; k < np; ++k) {
        if (t0[t] == k || t1[t] == k || t2[t] == k)
          continue;  // a plane that defines this vertex
        const Real s = n[k][0] * v[0] + n[k][1] * v[1] + n[k][2] * v[2] - nn[k];
        if (s >
            Real(0)) {  // v outside plane k; convert to a perpendicular distance for the tol test
          const Real invlen = (nn[k] > Real(0)) ? Real(1) / Kokkos::sqrt(nn[k]) : Real(0);
          maxd = Kokkos::max(maxd, s * invlen);
        }
      }
    }
    return maxd <= tol;
  }
 
  KOKKOS_INLINE_FUNCTION bool isSelfConsistentPartners(Real tol, int* outPartners, int maxPartners,
                                                       int& nPartners) const {
    nPartners = 0;
    bool consistent = true;
    for (int t = 0; t < nt; ++t) {
      if (!alive[t])
        continue;
      const Real v[3] = {vx[t], vy[t], vz[t]};
      for (int k = 0; k < np; ++k) {
        if (t0[t] == k || t1[t] == k || t2[t] == k)
          continue;  // a plane that defines this vertex
        const Real s = n[k][0] * v[0] + n[k][1] * v[1] + n[k][2] * v[2] - nn[k];
        if (s > Real(0)) {  // v outside plane k; perpendicular distance vs tol
          const Real invlen = (nn[k] > Real(0)) ? Real(1) / Kokkos::sqrt(nn[k]) : Real(0);
          if (s * invlen > tol) {
            consistent = false;
            const int part = pnbr[k];
            if (part >= 0) {  // box/SDF planes (pnbr<0) have no partner cell
              bool seen = false;
              for (int q = 0; q < nPartners; ++q)
                if (outPartners[q] == part) {
                  seen = true;
                  break;
                }
              if (!seen && nPartners < maxPartners)
                outPartners[nPartners++] = part;
            }
          }
        }
      }
    }
    return consistent;
  }
 
  KOKKOS_INLINE_FUNCTION int oppositePlane(int s, int x, int y) const {
    const int a = t0[s], b = t1[s], c = t2[s];
    return (a != x && a != y) ? a : (b != x && b != y) ? b : c;
  }
 
  KOKKOS_INLINE_FUNCTION void computePoke4(unsigned char* out) const {
    static_assert(TrackAdj, "computePoke4 requires TrackAdj");
    for (int t = 0; t < nt; ++t) {
      if (!alive[t]) {
        out[t * 4 + 0] = out[t * 4 + 1] = out[t * 4 + 2] = out[t * 4 + 3] = 0xff;
        continue;
      }
      const int vt[3] = {t0[t], t1[t], t2[t]};
      int eo[3];
      for (int e = 0; e < 3; ++e) {
        const int s = adj[t * 3 + e];
        eo[e] = (s >= 0 && s < nt && alive[s]) ? oppositePlane(s, vt[e], vt[(e + 1) % 3]) : 0xff;
        out[t * 4 + e] = (unsigned char)(eo[e] & 0xff);
      }
      const Real v[3] = {vx[t], vy[t], vz[t]};
      Real best = Real(3.4e38);
      int g = 0xff;
      for (int k = 0; k < np; ++k) {
        if (k == vt[0] || k == vt[1] || k == vt[2])
          continue;  // defining plane
        if (k == eo[0] || k == eo[1] || k == eo[2])
          continue;  // already an edge-opposite
        const Real slack = nn[k] - (n[k][0] * v[0] + n[k][1] * v[1] + n[k][2] * v[2]);
        if (slack >= Real(0) && slack < best) {
          best = slack;
          g = k;
        }
      }
      out[t * 4 + 3] = (unsigned char)(g & 0xff);
    }
  }
 
  KOKKOS_INLINE_FUNCTION bool isLocallyConvex(const unsigned char* poke4, Real tol) const {
    Real maxd = Real(0);
    for (int t = 0; t < nt; ++t) {
      if (!alive[t])
        continue;
      const Real v[3] = {vx[t], vy[t], vz[t]};
      for (int j = 0; j < 4; ++j) {
        const int k = (int)poke4[t * 4 + j];
        if (k == 0xff)
          continue;
        const Real sl = n[k][0] * v[0] + n[k][1] * v[1] + n[k][2] * v[2] - nn[k];
        if (sl > Real(0)) {
          const Real invlen = (nn[k] > Real(0)) ? Real(1) / Kokkos::sqrt(nn[k]) : Real(0);
          maxd = Kokkos::max(maxd, sl * invlen);
        }
      }
    }
    return maxd <= tol;
  }
 
  KOKKOS_INLINE_FUNCTION bool isLocallyConvexPartners(const unsigned char* poke4, Real tol,
                                                      int* outPartners, int maxPartners,
                                                      int& nPartners) const {
    nPartners = 0;
    bool consistent = true;
    for (int t = 0; t < nt; ++t) {
      if (!alive[t])
        continue;
      const Real v[3] = {vx[t], vy[t], vz[t]};
      for (int j = 0; j < 4; ++j) {
        const int k = (int)poke4[t * 4 + j];
        if (k == 0xff)
          continue;
        const Real sl = n[k][0] * v[0] + n[k][1] * v[1] + n[k][2] * v[2] - nn[k];
        if (sl > Real(0)) {
          const Real invlen = (nn[k] > Real(0)) ? Real(1) / Kokkos::sqrt(nn[k]) : Real(0);
          if (sl * invlen > tol) {
            consistent = false;
            const int part = pnbr[k];
            if (part >= 0) {
              bool seen = false;
              for (int q = 0; q < nPartners; ++q)
                if (outPartners[q] == part) {
                  seen = true;
                  break;
                }
              if (!seen && nPartners < maxPartners)
                outPartners[nPartners++] = part;
            }
          }
        }
      }
    }
    return consistent;
  }
 
  KOKKOS_INLINE_FUNCTION Real maxVertexRsq() const {
    Real m = 0;
    for (int t = 0; t < nt; ++t) {
      if (!alive[t])
        continue;
      const Real r = vx[t] * vx[t] + vy[t] * vy[t] + vz[t] * vz[t];
      if (r > m)
        m = r;
    }
    return m;
  }
 
  KOKKOS_INLINE_FUNCTION int findSharing(int self, int x, int y) const {
    for (int s = 0; s < nt; ++s) {
      if (s == self || !alive[s])
        continue;
      const int a = t0[s], b = t1[s], c = t2[s];
      const bool hasX = (a == x || b == x || c == x);
      const bool hasY = (a == y || b == y || c == y);
      if (hasX && hasY)
        return s;
    }
    return -1;
  }
 
  KOKKOS_INLINE_FUNCTION void rebuildAdjacency() {
    if constexpr (TrackAdj) {
      for (int t = 0; t < nt; ++t) {
        if (!alive[t])
          continue;
        const int vt[3] = {t0[t], t1[t], t2[t]};
        for (int e = 0; e < 3; ++e)
          adj[t * 3 + e] = findSharing(t, vt[e], vt[(e + 1) % 3]);
      }
    }
  }
 
  KOKKOS_INLINE_FUNCTION bool checkAdjacencyInvariant() const {
    if constexpr (!TrackAdj) {
      return true;
    } else {
      for (int t = 0; t < nt; ++t) {
        if (!alive[t])
          continue;
        const int vt[3] = {t0[t], t1[t], t2[t]};
        for (int e = 0; e < 3; ++e) {
          const int x = vt[e], y = vt[(e + 1) % 3];
          const int s = adj[t * 3 + e];
          if (s < 0 || s >= nt || !alive[s])
            return false;  // live target
          const int vs[3] = {t0[s], t1[s], t2[s]};
          const bool hasX = (vs[0] == x || vs[1] == x || vs[2] == x);
          const bool hasY = (vs[0] == y || vs[1] == y || vs[2] == y);
          if (!(hasX && hasY))
            return false;     // shares the edge
          bool back = false;  // symmetry
          for (int f = 0; f < 3; ++f)
            if (adj[s * 3 + f] == t)
              back = true;
          if (!back)
            return false;
        }
      }
      return true;
    }
  }
 
  KOKKOS_INLINE_FUNCTION int allocTri() {
    for (int s = 0; s < nt; ++s)
      if (!alive[s])
        return s;
    if (nt < MAXT)
      return nt++;
    overflow = true;
    return -1;
  }
 
  KOKKOS_INLINE_FUNCTION bool clip(const Real pdir[3], Real d, int nbr) {
    if (np >= MAXP) {
      overflow = true;
      return false;
    }
    const int pi = np;  // tentative index
    const Real l2 = pdir[0] * pdir[0] + pdir[1] * pdir[1] + pdir[2] * pdir[2];
    const Real a = (l2 > Real(0)) ? d / l2 : Real(0);  // foot scale: nf = a·pdir
    n[pi][0] = a * pdir[0];
    n[pi][1] = a * pdir[1];
    n[pi][2] = a * pdir[2];
    nn[pi] = n[pi][0] * n[pi][0] + n[pi][1] * n[pi][1] + n[pi][2] * n[pi][2];
    pnbr[pi] = nbr;
 
    // Mark triangles whose dual vertex is outside the new half-space (nf·v > nf·nf).
    bool kill[MAXT];
    bool any = false;
    for (int t = 0; t < nt; ++t) {
      kill[t] = false;
      if (!alive[t])
        continue;
      const Real s = n[pi][0] * vx[t] + n[pi][1] * vy[t] + n[pi][2] * vz[t] - nn[pi];
      if (s > 0) {
        kill[t] = true;
        any = true;
      }
    }
    if (!any)
      return false;  // candidate does not cut -> no-op, plane not committed
 
    np = pi + 1;  // commit the plane
 
    if constexpr (!TrackAdj) {
      // ---- Lean path (unchanged): collect the horizon by findSharing, then retriangulate. ----
      // For each killed triangle edge whose sharing triangle is alive, the new plane gets a
      // triangle (x, y, pi). Collect first, then mutate.
      unsigned char nA[MAXT], nB[MAXT];
      int nnew = 0;
      for (int t = 0; t < nt; ++t) {
        if (!kill[t])
          continue;
        const int vtx[3] = {t0[t], t1[t], t2[t]};
        for (int e = 0; e < 3; ++e) {
          const int x = vtx[e], y = vtx[(e + 1) % 3];
          const int other = findSharing(t, x, y);
          if (other >= 0 && !kill[other]) {
            if (nnew < MAXT) {
              nA[nnew] = (unsigned char)x;
              nB[nnew] = (unsigned char)y;
              ++nnew;
            } else {
              overflow = true;
            }
          }
        }
      }
      // Remove killed triangles.
      for (int t = 0; t < nt; ++t)
        if (kill[t])
          alive[t] = false;
      // Add the new fan around plane pi.
      for (int i = 0; i < nnew; ++i) {
        const int slot = allocTri();
        if (slot < 0)
          break;
        t0[slot] = nA[i];
        t1[slot] = nB[i];
        t2[slot] = (unsigned char)pi;
        alive[slot] = true;
        computeVertex(slot);
      }
    } else {
      // ---- TrackAdj path: walk the conflict-region boundary as ONE oriented cycle via `adj`
      // (O(conflict), no findSharing), then stitch the new fan's adjacency in O(1) per edge. The
      // conflict region (vertices outside the new half-space) is connected and its boundary is a
      // single cycle for a convex clip. ----
      int hA[MAXT], hB[MAXT], hSurv[MAXT],
          hBack[MAXT];  // horizon edge {A,B}, survivor across it, slot back
      int H = 0;
      // Start edge: a killed triangle with a surviving neighbour across some slot.
      int t0k = -1, e0 = -1;
      for (int t = 0; t < nt && t0k < 0; ++t) {
        if (!kill[t])
          continue;
        for (int e = 0; e < 3; ++e) {
          const int s = adj[t * 3 + e];
          if (s >= 0 && alive[s] && !kill[s]) {
            t0k = t;
            e0 = e;
            break;
          }
        }
      }
      if (t0k < 0) {
        overflow = true;
        return true;
      }  // no surviving boundary (cannot happen for a real cut)
      const int vstart[3] = {t0[t0k], t1[t0k], t2[t0k]};
      int A = vstart[e0], B = vstart[(e0 + 1) % 3];
      int tt = t0k, ee = e0;
      const int guardMax = 6 * nt + 12;
      int guard = 0;
      do {
        const int surv = adj[tt * 3 + ee];
        int back = 0;
        for (int f = 0; f < 3; ++f)
          if (adj[surv * 3 + f] == tt) {
            back = f;
            break;
          }
        if (H < MAXT) {
          hA[H] = A;
          hB[H] = B;
          hSurv[H] = surv;
          hBack[H] = back;
          ++H;
        } else {
          overflow = true;
          break;
        }
        // Pivot around vertex B, rotating through killed triangles, until we exit to the next
        // survivor.
        const int pivot = B;
        int rt = tt, arrival = ee, rotGuard = 0;
        while (true) {
          const int rv[3] = {t0[rt], t1[rt], t2[rt]};
          const int lv = (rv[0] == pivot)   ? 0
                         : (rv[1] == pivot) ? 1
                                            : 2;                      // local index of pivot in rt
          const int crossSlot = (arrival == lv) ? (lv + 2) % 3 : lv;  // the OTHER B-edge of rt
          const int rn = adj[rt * 3 + crossSlot];
          if (rn >= 0 && kill[rn]) {  // still inside the conflict region: keep rotating around B
            int f = (adj[rn * 3 + 0] == rt) ? 0 : (adj[rn * 3 + 1] == rt) ? 1 : 2;
            rt = rn;
            arrival = f;
            if (++rotGuard > guardMax) {
              overflow = true;
              break;
            }
            continue;
          }
          // rn survives -> edge crossSlot of rt is the next horizon edge, oriented out of the
          // pivot.
          const int u = rv[crossSlot], w = rv[(crossSlot + 1) % 3];
          A = pivot;
          B = (u == pivot) ? w : u;
          tt = rt;
          ee = crossSlot;
          break;
        }
        if (overflow)
          break;
        if (++guard > guardMax) {
          overflow = true;
          break;
        }
      } while (!(tt == t0k && ee == e0));
 
      // Remove killed triangles (frees their slots for the fan via allocTri).
      for (int t = 0; t < nt; ++t)
        if (kill[t])
          alive[t] = false;
      // Allocate the whole fan first (so τ_{i±1} cross-references resolve), then set topology +
      // adjacency.
      int tau[MAXT];
      bool ok = !overflow;
      for (int i = 0; i < H && ok; ++i) {
        tau[i] = allocTri();
        if (tau[i] < 0)
          ok = false;
        else
          alive[tau[i]] = true;  // reserve the slot so allocTri doesn't hand it out again
      }
      if (ok) {
        for (int i = 0; i < H; ++i) {
          const int s = tau[i];
          t0[s] = (unsigned char)hA[i];
          t1[s] = (unsigned char)hB[i];
          t2[s] = (unsigned char)pi;
          alive[s] = true;
          computeVertex(s);
          adj[s * 3 + 0] = hSurv[i];  // slot 0 = edge {A_i,B_i} -> survivor across the horizon
          adj[s * 3 + 1] = tau[(i + 1) % H];      // slot 1 = edge {B_i,P}   -> next fan triangle
          adj[s * 3 + 2] = tau[(i - 1 + H) % H];  // slot 2 = edge {P,A_i}   -> prev fan triangle
          adj[hSurv[i] * 3 + hBack[i]] = s;  // back-patch the survivor's slot to this fan triangle
        }
      }
    }
    return true;
  }
 
  KOKKOS_INLINE_FUNCTION bool empty() const {
    for (int t = 0; t < nt; ++t)
      if (alive[t])
        return false;
    return true;
  }
 
  KOKKOS_INLINE_FUNCTION int countFaces() const {
    int nf = 0;
    for (int k = 0; k < np; ++k) {
      bool used = false;
      for (int t = 0; t < nt && !used; ++t)
        if (alive[t] && (t0[t] == k || t1[t] == k || t2[t] == k))
          used = true;
      if (used)
        ++nf;
    }
    return nf;
  }
 
  static constexpr int MAXFV = 28;  // max vertices on one face polygon
 
  KOKKOS_INLINE_FUNCTION int faceOrdered(int k, Real fx[MAXFV], Real fy[MAXFV],
                                         Real fz[MAXFV]) const {
    int m = 0;
    for (int t = 0; t < nt; ++t) {
      if (!alive[t])
        continue;
      if (t0[t] != k && t1[t] != k && t2[t] != k)
        continue;
      if (m < MAXFV) {
        fx[m] = vx[t];
        fy[m] = vy[t];
        fz[m] = vz[t];
        ++m;
      }
    }
    if (m < 3)
      return m;
    const Real nx = n[k][0], ny = n[k][1], nz = n[k][2];
    const Real nlen = Kokkos::sqrt(nx * nx + ny * ny + nz * nz);
    if (nlen == Real(0))
      return 0;
    const Real un[3] = {nx / nlen, ny / nlen, nz / nlen};
    Real e1[3];
    if (Kokkos::fabs(un[0]) <= Kokkos::fabs(un[1]) && Kokkos::fabs(un[0]) <= Kokkos::fabs(un[2])) {
      e1[0] = 0;
      e1[1] = -un[2];
      e1[2] = un[1];
    } else if (Kokkos::fabs(un[1]) <= Kokkos::fabs(un[2])) {
      e1[0] = -un[2];
      e1[1] = 0;
      e1[2] = un[0];
    } else {
      e1[0] = -un[1];
      e1[1] = un[0];
      e1[2] = 0;
    }
    const Real e1l = Kokkos::sqrt(e1[0] * e1[0] + e1[1] * e1[1] + e1[2] * e1[2]);
    e1[0] /= e1l;
    e1[1] /= e1l;
    e1[2] /= e1l;
    const Real e2[3] = {un[1] * e1[2] - un[2] * e1[1], un[2] * e1[0] - un[0] * e1[2],
                        un[0] * e1[1] - un[1] * e1[0]};
    Real cx = 0, cy = 0, cz = 0;
    for (int i = 0; i < m; ++i) {
      cx += fx[i];
      cy += fy[i];
      cz += fz[i];
    }
    cx /= m;
    cy /= m;
    cz /= m;
    Real ang[MAXFV];
    for (int i = 0; i < m; ++i) {
      const Real dx = fx[i] - cx, dy = fy[i] - cy, dz = fz[i] - cz;
      const Real px = dx * e1[0] + dy * e1[1] + dz * e1[2];  // in-plane coords
      const Real py = dx * e2[0] + dy * e2[1] + dz * e2[2];
      // diamond pseudo-angle: monotonic in true angle, in [0,4), NO atan2 (transcendental)
      const Real s = Kokkos::fabs(px) + Kokkos::fabs(py);
      const Real t = (s > Real(0)) ? py / s : Real(0);  // [-1,1]
      ang[i] = (px < Real(0)) ? (Real(2) - t) : (py < Real(0) ? Real(4) + t : t);
    }
    for (int i = 1; i < m; ++i) {
      Real ka = ang[i], kx = fx[i], ky = fy[i], kz = fz[i];
      int j = i - 1;
      while (j >= 0 && ang[j] > ka) {
        ang[j + 1] = ang[j];
        fx[j + 1] = fx[j];
        fy[j + 1] = fy[j];
        fz[j + 1] = fz[j];
        --j;
      }
      ang[j + 1] = ka;
      fx[j + 1] = kx;
      fy[j + 1] = ky;
      fz[j + 1] = kz;
    }
    return m;
  }
 
  KOKKOS_INLINE_FUNCTION static void polyAreaVec(const Real fx[], const Real fy[], const Real fz[],
                                                 int m, Real A[3]) {
    Real ax = 0, ay = 0, az = 0;
    for (int i = 0; i < m; ++i) {
      const int j = (i + 1 == m) ? 0 : i + 1;
      ax += fy[i] * fz[j] - fz[i] * fy[j];
      ay += fz[i] * fx[j] - fx[i] * fz[j];
      az += fx[i] * fy[j] - fy[i] * fx[j];
    }
    A[0] = Real(0.5) * ax;
    A[1] = Real(0.5) * ay;
    A[2] = Real(0.5) * az;
  }
 
  static constexpr int MAXSV = 2 * MAXFV;  
 
  KOKKOS_INLINE_FUNCTION int sectionPolygon(const Real p0[3], const Real u[3], Real px[MAXSV],
                                            Real py[MAXSV], Real pz[MAXSV]) const {
    int m = 0;
    for (int i = 0; i < nt; ++i) {
      if (!alive[i])
        continue;
      const Real di = (vx[i] - p0[0]) * u[0] + (vy[i] - p0[1]) * u[1] + (vz[i] - p0[2]) * u[2];
      const unsigned char a0 = t0[i], a1 = t1[i], a2 = t2[i];
      for (int j = i + 1; j < nt; ++j) {
        if (!alive[j])
          continue;
        const unsigned char b0 = t0[j], b1 = t1[j], b2 = t2[j];
        const int common = (a0 == b0) + (a0 == b1) + (a0 == b2) + (a1 == b0) + (a1 == b1) +
                           (a1 == b2) + (a2 == b0) + (a2 == b1) + (a2 == b2);
        if (common != 2)
          continue;  // not an edge (a primal edge = two triangles sharing exactly two planes)
        const Real dj = (vx[j] - p0[0]) * u[0] + (vy[j] - p0[1]) * u[1] + (vz[j] - p0[2]) * u[2];
        if ((di > Real(0)) == (dj > Real(0)))
          continue;  // both endpoints on the same side -> edge doesn't cross the plane
        const Real denom = di - dj;
        if (denom == Real(0))
          continue;
        const Real tt = di / denom;
        if (m >= MAXSV)
          return -1;
        px[m] = vx[i] + tt * (vx[j] - vx[i]);
        py[m] = vy[i] + tt * (vy[j] - vy[i]);
        pz[m] = vz[i] + tt * (vz[j] - vz[i]);
        ++m;
      }
    }
    if (m < 3)
      return m;
    const Real ulen = Kokkos::sqrt(u[0] * u[0] + u[1] * u[1] + u[2] * u[2]);
    if (ulen == Real(0))
      return 0;
    const Real un[3] = {u[0] / ulen, u[1] / ulen, u[2] / ulen};
    Real e1[3];
    if (Kokkos::fabs(un[0]) <= Kokkos::fabs(un[1]) && Kokkos::fabs(un[0]) <= Kokkos::fabs(un[2])) {
      e1[0] = 0; e1[1] = -un[2]; e1[2] = un[1];
    } else if (Kokkos::fabs(un[1]) <= Kokkos::fabs(un[2])) {
      e1[0] = -un[2]; e1[1] = 0; e1[2] = un[0];
    } else {
      e1[0] = -un[1]; e1[1] = un[0]; e1[2] = 0;
    }
    const Real e1l = Kokkos::sqrt(e1[0] * e1[0] + e1[1] * e1[1] + e1[2] * e1[2]);
    e1[0] /= e1l; e1[1] /= e1l; e1[2] /= e1l;
    const Real e2[3] = {un[1] * e1[2] - un[2] * e1[1], un[2] * e1[0] - un[0] * e1[2],
                        un[0] * e1[1] - un[1] * e1[0]};
    Real cx = 0, cy = 0, cz = 0;
    for (int i = 0; i < m; ++i) { cx += px[i]; cy += py[i]; cz += pz[i]; }
    cx /= m; cy /= m; cz /= m;
    Real ang[MAXSV];
    for (int i = 0; i < m; ++i) {
      const Real dx = px[i] - cx, dy = py[i] - cy, dz = pz[i] - cz;
      const Real qx = dx * e1[0] + dy * e1[1] + dz * e1[2];
      const Real qy = dx * e2[0] + dy * e2[1] + dz * e2[2];
      const Real s = Kokkos::fabs(qx) + Kokkos::fabs(qy);
      const Real tt = (s > Real(0)) ? qy / s : Real(0);
      ang[i] = (qx < Real(0)) ? (Real(2) - tt) : (qy < Real(0) ? Real(4) + tt : tt);
    }
    for (int i = 1; i < m; ++i) {
      Real ka = ang[i], kx = px[i], ky = py[i], kz = pz[i];
      int j = i - 1;
      while (j >= 0 && ang[j] > ka) {
        ang[j + 1] = ang[j];
        px[j + 1] = px[j]; py[j + 1] = py[j]; pz[j + 1] = pz[j];
        --j;
      }
      ang[j + 1] = ka; px[j + 1] = kx; py[j + 1] = ky; pz[j + 1] = kz;
    }
    return m;
  }
 
  KOKKOS_INLINE_FUNCTION Real volume() const {
    Real vol = 0;
    Real fx[MAXFV], fy[MAXFV], fz[MAXFV];
    for (int k = 0; k < np; ++k) {
      const int m = faceOrdered(k, fx, fy, fz);
      if (m < 3)
        continue;
      Real A[3];
      polyAreaVec(fx, fy, fz, m, A);
      const Real area = Kokkos::sqrt(A[0] * A[0] + A[1] * A[1] + A[2] * A[2]);
      const Real support = Kokkos::sqrt(nn[k]);  // |n_k| = seed->plane distance
      vol += support * area;
    }
    return vol * (Real(1) / Real(3));
  }
 
  // ---- Vertex-local, SORT-FREE geometry (Ray/Sokolov/Lefebvre/Lévy TOG 2018; geogram ConvexCell)
  // ---- Plane stored as the (non-unit) foot-point normal n with interior {x : n·x ≤ n·n}; then x=n
  // is the foot of the perpendicular from the origin onto the plane (the facet point), and |n| is
  // the origin→ plane distance. The cell now stores n directly (member `n`, with `nn = n·n`), so
  // this is just a read. Volume by the divergence theorem, coning every boundary flag (facet n_i,
  // edge foot f, vertex v) to the origin: each tetra (0,n_i,f,v) is LOCAL to one vertex + its 3
  // planes, and the signed sum is exact even when feet fall outside their faces. No vertex
  // ordering, no adjacency — a pure per-vertex scatter.
 
  KOKKOS_INLINE_FUNCTION void planeN(int k, Real out[3]) const {
    out[0] = n[k][0];
    out[1] = n[k][1];
    out[2] = n[k][2];
  }
  KOKKOS_INLINE_FUNCTION static void xprod(const Real a[3], const Real b[3], Real o[3]) {
    o[0] = a[1] * b[2] - a[2] * b[1];
    o[1] = a[2] * b[0] - a[0] * b[2];
    o[2] = a[0] * b[1] - a[1] * b[0];
  }
  KOKKOS_INLINE_FUNCTION static Real dot3(const Real a[3], const Real b[3]) {
    return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
  }
  KOKKOS_INLINE_FUNCTION static Real det3(const Real a[3], const Real b[3], const Real c[3]) {
    Real bc[3];
    xprod(b, c, bc);
    return dot3(a, bc);
  }
  KOKKOS_INLINE_FUNCTION static void edgeFoot(const Real v[3], const Real ck[3], Real f[3]) {
    const Real c2 = dot3(ck, ck), s = (c2 > Real(0)) ? dot3(v, ck) / c2 : Real(0);
    f[0] = v[0] - s * ck[0];
    f[1] = v[1] - s * ck[1];
    f[2] = v[2] - s * ck[2];
  }
 
  KOKKOS_INLINE_FUNCTION Real volumePerVertex() const {
    Real vol = 0;
    for (int t = 0; t < nt; ++t) {
      if (!alive[t])
        continue;
      Real n1[3], n2[3], n3[3];
      planeN(t0[t], n1);
      planeN(t1[t], n2);
      planeN(t2[t], n3);
      Real c1[3], c2[3], c3[3];
      xprod(n2, n3, c1);
      xprod(n3, n1, c2);
      xprod(n1, n2, c3);
      Real D = dot3(n1, c1);
      if (D < Real(0))  // canonical order D>0: swap (n2,n3) and (c2,c3); v is unchanged
        for (int a = 0; a < 3; ++a) {
          Real tmp = n2[a];
          n2[a] = n3[a];
          n3[a] = tmp;
          tmp = c2[a];
          c2[a] = c3[a];
          c3[a] = tmp;
        }
      const Real v[3] = {vx[t], vy[t], vz[t]};
      Real f12[3], f23[3], f31[3];
      edgeFoot(v, c3, f12);
      edgeFoot(v, c1, f23);
      edgeFoot(v, c2, f31);
      Real e[3], d = 0;
      e[0] = n1[0] - n2[0];
      e[1] = n1[1] - n2[1];
      e[2] = n1[2] - n2[2];
      d += det3(e, f12, v);
      e[0] = n2[0] - n3[0];
      e[1] = n2[1] - n3[1];
      e[2] = n2[2] - n3[2];
      d += det3(e, f23, v);
      e[0] = n3[0] - n1[0];
      e[1] = n3[1] - n1[1];
      e[2] = n3[2] - n1[2];
      d += det3(e, f31, v);
      vol += d;
    }
    return vol * (Real(1) / Real(6));
  }
 
  KOKKOS_INLINE_FUNCTION void facetAreasPerVertex(Real* area) const {
    for (int t = 0; t < nt; ++t) {
      if (!alive[t])
        continue;
      int k1 = t0[t], k2 = t1[t], k3 = t2[t];
      Real n1[3], n2[3], n3[3];
      planeN(k1, n1);
      planeN(k2, n2);
      planeN(k3, n3);
      Real c1[3], c2[3], c3[3];
      xprod(n2, n3, c1);
      xprod(n3, n1, c2);
      xprod(n1, n2, c3);
      Real D = dot3(n1, c1);
      if (D < Real(0)) {  // canonical: swap planes 2,3 (normals, cross, AND indices)
        for (int a = 0; a < 3; ++a) {
          Real tm = n2[a];
          n2[a] = n3[a];
          n3[a] = tm;
          tm = c2[a];
          c2[a] = c3[a];
          c3[a] = tm;
        }
        int tk = k2;
        k2 = k3;
        k3 = tk;
      }
      const Real v[3] = {vx[t], vy[t], vz[t]};
      Real f12[3], f23[3], f31[3];
      edgeFoot(v, c3, f12);
      edgeFoot(v, c1, f23);
      edgeFoot(v, c2, f31);
      Real g[3];
      g[0] = f12[0] - f31[0];
      g[1] = f12[1] - f31[1];
      g[2] = f12[2] - f31[2];
      area[k1] += det3(n1, g, v) / (Real(2) * Kokkos::sqrt(dot3(n1, n1)));
      g[0] = f23[0] - f12[0];
      g[1] = f23[1] - f12[1];
      g[2] = f23[2] - f12[2];
      area[k2] += det3(n2, g, v) / (Real(2) * Kokkos::sqrt(dot3(n2, n2)));
      g[0] = f31[0] - f23[0];
      g[1] = f31[1] - f23[1];
      g[2] = f31[2] - f23[2];
      area[k3] += det3(n3, g, v) / (Real(2) * Kokkos::sqrt(dot3(n3, n3)));
    }
  }
 
  KOKKOS_INLINE_FUNCTION static void scatterFacetMoment(int ki, const Real ni[3], const Real ff[3],
                                                        const Real fl[3], const Real v[3],
                                                        Real* area, Real* mx, Real* my, Real* mz) {
    const Real inv2n = Real(0.5) / Kokkos::sqrt(dot3(ni, ni));
    const Real sa = det3(ni, ff, v) * inv2n;   // triangle (ni, ff, v)
    const Real sb = -det3(ni, fl, v) * inv2n;  // triangle (ni, v, fl)
    area[ki] += sa + sb;
    mx[ki] += sa * (ni[0] + ff[0] + v[0]) / Real(3) + sb * (ni[0] + v[0] + fl[0]) / Real(3);
    my[ki] += sa * (ni[1] + ff[1] + v[1]) / Real(3) + sb * (ni[1] + v[1] + fl[1]) / Real(3);
    mz[ki] += sa * (ni[2] + ff[2] + v[2]) / Real(3) + sb * (ni[2] + v[2] + fl[2]) / Real(3);
  }
 
  KOKKOS_INLINE_FUNCTION void facetMomentsPerVertex(Real* area, Real* mx, Real* my,
                                                    Real* mz) const {
    for (int t = 0; t < nt; ++t) {
      if (!alive[t])
        continue;
      int k1 = t0[t], k2 = t1[t], k3 = t2[t];
      Real n1[3], n2[3], n3[3];
      planeN(k1, n1);
      planeN(k2, n2);
      planeN(k3, n3);
      Real c1[3], c2[3], c3[3];
      xprod(n2, n3, c1);
      xprod(n3, n1, c2);
      xprod(n1, n2, c3);
      if (dot3(n1, c1) < Real(0)) {  // canonical: swap planes 2,3 (normals, cross, indices)
        for (int a = 0; a < 3; ++a) {
          Real tm = n2[a];
          n2[a] = n3[a];
          n3[a] = tm;
          tm = c2[a];
          c2[a] = c3[a];
          c3[a] = tm;
        }
        int tk = k2;
        k2 = k3;
        k3 = tk;
      }
      const Real v[3] = {vx[t], vy[t], vz[t]};
      Real f12[3], f23[3], f31[3];
      edgeFoot(v, c3, f12);
      edgeFoot(v, c1, f23);
      edgeFoot(v, c2, f31);
      scatterFacetMoment(k1, n1, f12, f31, v, area, mx, my, mz);
      scatterFacetMoment(k2, n2, f23, f12, v, area, mx, my, mz);
      scatterFacetMoment(k3, n3, f31, f23, v, area, mx, my, mz);
    }
  }
 
  KOKKOS_INLINE_FUNCTION static void scatterGradRaw(int ki, const Real ni[3], const Real ff[3],
                                                    const Real fl[3], const Real v[3], Real* dgx,
                                                    Real* dgy, Real* dgz) {
    const Real sa = det3(ni, ff, v);   // raw signed area of triangle (ni, ff, v)
    const Real sb = -det3(ni, fl, v);  // raw signed area of triangle (ni, v, fl)
    const Real da = sa + sb;
    const Real momx = sa * (ni[0] + ff[0] + v[0]) / Real(3) + sb * (ni[0] + v[0] + fl[0]) / Real(3);
    const Real momy = sa * (ni[1] + ff[1] + v[1]) / Real(3) + sb * (ni[1] + v[1] + fl[1]) / Real(3);
    const Real momz = sa * (ni[2] + ff[2] + v[2]) / Real(3) + sb * (ni[2] + v[2] + fl[2]) / Real(3);
    dgx[ki] += Real(2) * da * ni[0] - momx;  // 2·da·n_k − mom  (numerator piece)
    dgy[ki] += Real(2) * da * ni[1] - momy;
    dgz[ki] += Real(2) * da * ni[2] - momz;
  }
 
  KOKKOS_INLINE_FUNCTION void geomVolumeGrad(Real& vol, Real* dgx, Real* dgy, Real* dgz) const {
    Real V = 0;
    for (int t = 0; t < nt; ++t) {
      if (!alive[t])
        continue;
      int k1 = t0[t], k2 = t1[t], k3 = t2[t];
      Real n1[3], n2[3], n3[3];
      planeN(k1, n1);
      planeN(k2, n2);
      planeN(k3, n3);
      Real c1[3], c2[3], c3[3];
      xprod(n2, n3, c1);
      xprod(n3, n1, c2);
      xprod(n1, n2, c3);
      if (dot3(n1, c1) < Real(0)) {  // canonical D>0: swap planes 2,3 (normals, cross, indices)
        for (int a = 0; a < 3; ++a) {
          Real tm = n2[a];
          n2[a] = n3[a];
          n3[a] = tm;
          tm = c2[a];
          c2[a] = c3[a];
          c3[a] = tm;
        }
        int tk = k2;
        k2 = k3;
        k3 = tk;
      }
      const Real v[3] = {vx[t], vy[t], vz[t]};
      Real f12[3], f23[3], f31[3];
      edgeFoot(v, c3, f12);
      edgeFoot(v, c1, f23);
      edgeFoot(v, c2, f31);
      Real e[3];
      e[0] = n1[0] - n2[0];
      e[1] = n1[1] - n2[1];
      e[2] = n1[2] - n2[2];
      V += det3(e, f12, v);
      e[0] = n2[0] - n3[0];
      e[1] = n2[1] - n3[1];
      e[2] = n2[2] - n3[2];
      V += det3(e, f23, v);
      e[0] = n3[0] - n1[0];
      e[1] = n3[1] - n1[1];
      e[2] = n3[2] - n1[2];
      V += det3(e, f31, v);
      scatterGradRaw(k1, n1, f12, f31, v, dgx, dgy, dgz);
      scatterGradRaw(k2, n2, f23, f12, v, dgx, dgy, dgz);
      scatterGradRaw(k3, n3, f31, f23, v, dgx, dgy, dgz);
    }
    for (int k = 0; k < np; ++k) {  // per-facet fold, sqrt-free: numerator/(2·nn)
      const Real inv = (nn[k] > Real(0)) ? Real(1) / (Real(2) * nn[k]) : Real(0);
      dgx[k] *= inv;
      dgy[k] *= inv;
      dgz[k] *= inv;
    }
    vol = V * (Real(1) / Real(6));
  }
 
  KOKKOS_INLINE_FUNCTION static void scatterAreaMomentRaw(int ki, const Real ni[3],
                                                          const Real ff[3], const Real fl[3],
                                                          const Real v[3], Real* sa_acc, Real* smx,
                                                          Real* smy, Real* smz) {
    const Real sa = det3(ni, ff, v);   // raw signed area of triangle (ni, ff, v)
    const Real sb = -det3(ni, fl, v);  // raw signed area of triangle (ni, v, fl)
    sa_acc[ki] += sa + sb;
    smx[ki] += sa * (ni[0] + ff[0] + v[0]) / Real(3) + sb * (ni[0] + v[0] + fl[0]) / Real(3);
    smy[ki] += sa * (ni[1] + ff[1] + v[1]) / Real(3) + sb * (ni[1] + v[1] + fl[1]) / Real(3);
    smz[ki] += sa * (ni[2] + ff[2] + v[2]) / Real(3) + sb * (ni[2] + v[2] + fl[2]) / Real(3);
  }
 
  KOKKOS_INLINE_FUNCTION void geomVolumeArea(Real& vol, Real* avx, Real* avy, Real* avz, Real* dgx,
                                             Real* dgy, Real* dgz) const {
    Real V = 0;
    for (int t = 0; t < nt; ++t) {
      if (!alive[t])
        continue;
      int k1 = t0[t], k2 = t1[t], k3 = t2[t];
      Real n1[3], n2[3], n3[3];
      planeN(k1, n1);
      planeN(k2, n2);
      planeN(k3, n3);
      Real c1[3], c2[3], c3[3];
      xprod(n2, n3, c1);
      xprod(n3, n1, c2);
      xprod(n1, n2, c3);
      if (dot3(n1, c1) < Real(0)) {  // canonical D>0: swap planes 2,3 (normals, cross, indices)
        for (int a = 0; a < 3; ++a) {
          Real tm = n2[a];
          n2[a] = n3[a];
          n3[a] = tm;
          tm = c2[a];
          c2[a] = c3[a];
          c3[a] = tm;
        }
        int tk = k2;
        k2 = k3;
        k3 = tk;
      }
      const Real v[3] = {vx[t], vy[t], vz[t]};
      Real f12[3], f23[3], f31[3];
      edgeFoot(v, c3, f12);
      edgeFoot(v, c1, f23);
      edgeFoot(v, c2, f31);
      Real e[3];
      e[0] = n1[0] - n2[0];
      e[1] = n1[1] - n2[1];
      e[2] = n1[2] - n2[2];
      V += det3(e, f12, v);
      e[0] = n2[0] - n3[0];
      e[1] = n2[1] - n3[1];
      e[2] = n2[2] - n3[2];
      V += det3(e, f23, v);
      e[0] = n3[0] - n1[0];
      e[1] = n3[1] - n1[1];
      e[2] = n3[2] - n1[2];
      V += det3(e, f31, v);
      scatterAreaMomentRaw(k1, n1, f12, f31, v, avx, dgx, dgy,
                           dgz);  // avx <- S_a, (dgx,dgy,dgz) <- S_m
      scatterAreaMomentRaw(k2, n2, f23, f12, v, avx, dgx, dgy, dgz);
      scatterAreaMomentRaw(k3, n3, f31, f23, v, avx, dgx, dgy, dgz);
    }
    for (int k = 0; k < np; ++k) {  // per-facet fold, sqrt-free
      const Real inv = (nn[k] > Real(0)) ? Real(1) / (Real(2) * nn[k]) : Real(0);
      const Real sa = avx[k];  // raw area S_a (accumulated here during the scatter)
      const Real smx = dgx[k], smy = dgy[k], smz = dgz[k];
      dgx[k] = (Real(2) * sa * n[k][0] - smx) * inv;  // dV/dn_k
      dgy[k] = (Real(2) * sa * n[k][1] - smy) * inv;
      dgz[k] = (Real(2) * sa * n[k][2] - smz) * inv;
      avx[k] = sa * n[k][0] * inv;  // outward area-vector A_k·n_k/|n_k|
      avy[k] = sa * n[k][1] * inv;
      avz[k] = sa * n[k][2] * inv;
    }
    vol = V * (Real(1) / Real(6));
  }
 
  KOKKOS_INLINE_FUNCTION static void edgeCrossAdj(int j, int p, int q, const Real Nn[3][3],
                                                  const Real w[3], Real out[3]) {
    if (j == p) {
      xprod(Nn[q], w, out);
    } else if (j == q) {
      Real tm[3];
      xprod(Nn[p], w, tm);
      out[0] = -tm[0];
      out[1] = -tm[1];
      out[2] = -tm[2];
    } else {
      out[0] = out[1] = out[2] = Real(0);
    }
  }
 
  KOKKOS_INLINE_FUNCTION void geomVolumeAreaGrad(int t, int pl[3], Real contrib[3],
                                                 Real grad[3][3][3]) const {
    // planes + canonical swap (matches facetAreasPerVertex), reusing the cached vertex — nothing
    // recomputed. Every loop is VOR_UNROLL'd so the small local arrays scalarize into registers on
    // GPU (no local-memory spill from runtime indices); on CPU the macro is empty.
    int kk[3] = {t0[t], t1[t], t2[t]};
    Real Nn[3][3];
    VOR_UNROLL for (int a = 0; a < 3; ++a) {
      Nn[a][0] = n[kk[a]][0];
      Nn[a][1] = n[kk[a]][1];
      Nn[a][2] = n[kk[a]][2];
    }
    {
      Real c01[3];
      xprod(Nn[1], Nn[2], c01);
      if (dot3(Nn[0], c01) < Real(0)) {
        VOR_UNROLL for (int d = 0; d < 3; ++d) {
          Real tm = Nn[1][d];
          Nn[1][d] = Nn[2][d];
          Nn[2][d] = tm;
        }
        int tk = kk[1];
        kk[1] = kk[2];
        kk[2] = tk;
      }
    }
    pl[0] = kk[0];
    pl[1] = kk[1];
    pl[2] = kk[2];
    const Real v[3] = {vx[t], vy[t], vz[t]};
    Real nn3[3], u[3][3];  // u[a] = 2 n_a − v
    VOR_UNROLL for (int a = 0; a < 3; ++a) {
      nn3[a] = dot3(Nn[a], Nn[a]);
      VOR_UNROLL for (int d = 0; d < 3; ++d) u[a][d] = Real(2) * Nn[a][d] - v[d];
    }
    // edge directions e[m] = N[eP]×N[eQ] for edges (0,1),(1,2),(2,0); cofactor of plane j is c_j =
    // e[(j+1)%3]
    const int eP[3] = {0, 1, 2}, eQ[3] = {1, 2, 0};
    Real e[3][3], ecc[3], es[3], F[3][3];
    VOR_UNROLL for (int m = 0; m < 3; ++m) {
      xprod(Nn[eP[m]], Nn[eQ[m]], e[m]);
      ecc[m] = dot3(e[m], e[m]);
      const Real evc = dot3(v, e[m]);
      es[m] = (ecc[m] > Real(0)) ? evc / ecc[m] : Real(0);
      VOR_UNROLL for (int d = 0; d < 3; ++d) F[m][d] =
          v[d] - es[m] * e[m][d];  // foot of v on edge line m
    }
    const Real D = dot3(Nn[0], e[1]);  // e[1] = N1×N2 = cofactor c_0;  D = det(n0,n1,n2) > 0
    Real g[3][3], den[3];
    VOR_UNROLL for (int i = 0; i < 3; ++i) {  // facet i: g_i = F[i] − F[(i+2)%3]
      const int mb = (i + 2) % 3;
      VOR_UNROLL for (int d = 0; d < 3; ++d) g[i][d] = F[i][d] - F[mb][d];
      const Real Ni = det3(Nn[i], g[i], v);
      den[i] = Real(2) * Kokkos::sqrt(nn3[i]);
      contrib[i] = (den[i] > Real(0)) ? Ni / den[i] : Real(0);
    }
    VOR_UNROLL for (int i = 0; i < 3; ++i) {
      Real wi[3];
      xprod(v, Nn[i], wi);  // w_i = v×n_i  (∂det3(n_i,·,v)/∂g)
      Real nixg[3];
      xprod(Nn[i], g[i], nixg);  // n_i×g_i      (∂det3(n_i,g_i,·)/∂v)
      const int mb = (i + 2) % 3;
      VOR_UNROLL for (int j = 0; j < 3; ++j) {
        const Real* cj = e[(j + 1) % 3];  // cofactor of plane j
        Real gN[3] = {0, 0, 0};
        if (i == j) {
          Real gxv[3];
          xprod(g[i], v, gxv);
          VOR_UNROLL for (int d = 0; d < 3; ++d) gN[d] += gxv[d];
        }  // ∂/∂n_i
        const Real alpha = dot3(nixg, cj) / D;  // ∂/∂v (rank-1)
        VOR_UNROLL for (int d = 0; d < 3; ++d) gN[d] += alpha * u[j][d];
        // ∂/∂g = (∂F[i]/∂n_j − ∂F[(i+2)%3]/∂n_j)^T w_i  — the two foot adjoints
        VOR_UNROLL for (int s = 0; s < 2; ++s) {
          const int m = (s == 0) ? i : mb;
          const Real sgn = (s == 0) ? Real(1) : Real(-1);
          const int p = eP[m], q = eQ[m];
          const Real cjw = dot3(cj, wi);    // for vAdj
          const Real cje = dot3(cj, e[m]);  // for dS
          const Real emw = dot3(e[m], wi);
          Real eaW[3], eaV[3], eaE[3];
          edgeCrossAdj(j, p, q, Nn, wi, eaW);
          edgeCrossAdj(j, p, q, Nn, v, eaV);
          edgeCrossAdj(j, p, q, Nn, e[m], eaE);
          const Real inv_ecc = (ecc[m] > Real(0)) ? Real(1) / ecc[m] : Real(0);
          Real dS[3];  // ∂s_m/∂n_j
          VOR_UNROLL for (int d = 0; d < 3; ++d) dS[d] =
              ((cje / D) * u[j][d] + eaV[d]) * inv_ecc - (es[m] * inv_ecc) * Real(2) * eaE[d];
          // (∂F[m]/∂n_j)^T w_i = (c_j·w_i/D) u_j − (e_m·w_i) dS − s_m (∂e_m/∂n_j)^T w_i
          VOR_UNROLL for (int d = 0; d < 3; ++d) gN[d] +=
              sgn * ((cjw / D) * u[j][d] - emw * dS[d] - es[m] * eaW[d]);
        }
        const Real invden = (den[i] > Real(0)) ? Real(1) / den[i] : Real(0);
        VOR_UNROLL for (int d = 0; d < 3; ++d) grad[i][j][d] = gN[d] * invden;
        if (i == j) {  // 1/(2|n_i|) factor: −A_i n_i/nn_i
          const Real f = (nn3[i] > Real(0)) ? contrib[i] / nn3[i] : Real(0);
          VOR_UNROLL for (int d = 0; d < 3; ++d) grad[i][j][d] -= f * Nn[i][d];
        }
      }
    }
  }
 
  KOKKOS_INLINE_FUNCTION void geomVolumeAreaGradAD(int t, int pl[3], Real contrib[3],
                                                   Real grad[3][3][3]) const {
    using D = detail::Dual<Real, 9>;
    int k1 = t0[t], k2 = t1[t], k3 = t2[t];
    Real rn1[3] = {n[k1][0], n[k1][1], n[k1][2]};
    Real rn2[3] = {n[k2][0], n[k2][1], n[k2][2]};
    Real rn3[3] = {n[k3][0], n[k3][1], n[k3][2]};
    {  // canonical D>0 swap on the real values (same as facetAreasPerVertex) so contributions land
       // right
      Real cc[3];
      xprod(rn2, rn3, cc);
      if (dot3(rn1, cc) < Real(0)) {
        for (int a = 0; a < 3; ++a) {
          Real tm = rn2[a];
          rn2[a] = rn3[a];
          rn3[a] = tm;
        }
        int tk = k2;
        k2 = k3;
        k3 = tk;
      }
    }
    pl[0] = k1;
    pl[1] = k2;
    pl[2] = k3;
    D n1[3], n2[3], n3[3];
    for (int a = 0; a < 3; ++a) {  // seed: n1 <- slots 0..2, n2 <- 3..5, n3 <- 6..8
      n1[a] = detail::dseed<Real, 9>(rn1[a], a);
      n2[a] = detail::dseed<Real, 9>(rn2[a], 3 + a);
      n3[a] = detail::dseed<Real, 9>(rn3[a], 6 + a);
    }
    const D nn1 = detail::ddot<Real, 9>(n1, n1), nn2 = detail::ddot<Real, 9>(n2, n2),
            nn3 = detail::ddot<Real, 9>(n3, n3);
    D c1[3], c2[3], c3[3];
    detail::dcross<Real, 9>(n2, n3, c1);
    detail::dcross<Real, 9>(n3, n1, c2);
    detail::dcross<Real, 9>(n1, n2, c3);
    const D Dd = detail::ddot<Real, 9>(n1, c1);  // det
    D v[3];
    for (int a = 0; a < 3; ++a)
      v[a] = (nn1 * c1[a] + nn2 * c2[a] + nn3 * c3[a]) / Dd;  // Cramer vertex
    D f12[3], f23[3], f31[3];
    detail::dedgeFoot<Real, 9>(v, c3, f12);
    detail::dedgeFoot<Real, 9>(v, c1, f23);
    detail::dedgeFoot<Real, 9>(v, c2, f31);
    const D two = detail::dnum<Real, 9>(Real(2));
    D g[3], aa[3];
    for (int a = 0; a < 3; ++a)
      g[a] = f12[a] - f31[a];
    aa[0] = detail::ddet3<Real, 9>(n1, g, v) / (two * detail::dsqrt<Real, 9>(nn1));
    for (int a = 0; a < 3; ++a)
      g[a] = f23[a] - f12[a];
    aa[1] = detail::ddet3<Real, 9>(n2, g, v) / (two * detail::dsqrt<Real, 9>(nn2));
    for (int a = 0; a < 3; ++a)
      g[a] = f31[a] - f23[a];
    aa[2] = detail::ddet3<Real, 9>(n3, g, v) / (two * detail::dsqrt<Real, 9>(nn3));
    for (int i = 0; i < 3; ++i) {
      contrib[i] = aa[i].v;
      for (int j = 0; j < 3; ++j)
        for (int c = 0; c < 3; ++c)
          grad[i][j][c] = aa[i].d[3 * j + c];
    }
  }
 
  KOKKOS_INLINE_FUNCTION void geometryPerVertex(Real& vol, Real* area, Real* mx, Real* my,
                                                Real* mz) const {
    Real V = 0;
    for (int t = 0; t < nt; ++t) {
      if (!alive[t])
        continue;
      int k1 = t0[t], k2 = t1[t], k3 = t2[t];
      Real n1[3], n2[3], n3[3];
      planeN(k1, n1);
      planeN(k2, n2);
      planeN(k3, n3);
      Real c1[3], c2[3], c3[3];
      xprod(n2, n3, c1);
      xprod(n3, n1, c2);
      xprod(n1, n2, c3);
      if (dot3(n1, c1) < Real(0)) {  // canonical D>0: swap planes 2,3 (normals, cross, indices)
        for (int a = 0; a < 3; ++a) {
          Real tm = n2[a];
          n2[a] = n3[a];
          n3[a] = tm;
          tm = c2[a];
          c2[a] = c3[a];
          c3[a] = tm;
        }
        int tk = k2;
        k2 = k3;
        k3 = tk;
      }
      const Real v[3] = {vx[t], vy[t], vz[t]};
      Real f12[3], f23[3], f31[3];
      edgeFoot(v, c3, f12);
      edgeFoot(v, c1, f23);
      edgeFoot(v, c2, f31);
      Real e[3];
      e[0] = n1[0] - n2[0];
      e[1] = n1[1] - n2[1];
      e[2] = n1[2] - n2[2];
      V += det3(e, f12, v);
      e[0] = n2[0] - n3[0];
      e[1] = n2[1] - n3[1];
      e[2] = n2[2] - n3[2];
      V += det3(e, f23, v);
      e[0] = n3[0] - n1[0];
      e[1] = n3[1] - n1[1];
      e[2] = n3[2] - n1[2];
      V += det3(e, f31, v);
      scatterFacetMoment(k1, n1, f12, f31, v, area, mx, my, mz);
      scatterFacetMoment(k2, n2, f23, f12, v, area, mx, my, mz);
      scatterFacetMoment(k3, n3, f31, f23, v, area, mx, my, mz);
    }
    vol = V * (Real(1) / Real(6));
  }
 
  KOKKOS_INLINE_FUNCTION bool facetGeometry(int k, Real areaVec[3], Real dv[3],
                                            Real conn[3]) const {
    Real fx[MAXFV], fy[MAXFV], fz[MAXFV];
    const int m = faceOrdered(k, fx, fy, fz);
    if (m < 3)
      return false;
    Real A[3];
    polyAreaVec(fx, fy, fz, m, A);
    const Real r[3] = {Real(2) * n[k][0], Real(2) * n[k][1],
                       Real(2) * n[k][2]};                    // connector r = 2·foot
    if (A[0] * r[0] + A[1] * r[1] + A[2] * r[2] < Real(0)) {  // orient outward (toward neighbour)
      A[0] = -A[0];
      A[1] = -A[1];
      A[2] = -A[2];
    }
    const Real area = Kokkos::sqrt(A[0] * A[0] + A[1] * A[1] + A[2] * A[2]);
    // area-weighted centroid (fan-triangulate from v0)
    Real cx = 0, cy = 0, cz = 0, asum = 0;
    for (int i = 1; i + 1 < m; ++i) {
      const Real ux = fx[i] - fx[0], uy = fy[i] - fy[0], uz = fz[i] - fz[0];
      const Real wx = fx[i + 1] - fx[0], wy = fy[i + 1] - fy[0], wz = fz[i + 1] - fz[0];
      const Real tx = uy * wz - uz * wy, ty = uz * wx - ux * wz, tz = ux * wy - uy * wx;
      const Real at = Real(0.5) * Kokkos::sqrt(tx * tx + ty * ty + tz * tz);
      cx += at * (fx[0] + fx[i] + fx[i + 1]);
      cy += at * (fy[0] + fy[i] + fy[i + 1]);
      cz += at * (fz[0] + fz[i] + fz[i + 1]);
      asum += at;
    }
    const Real inv = asum > Real(0) ? Real(1) / (Real(3) * asum) : Real(0);
    cx *= inv;
    cy *= inv;
    cz *= inv;
    const Real s = area / Kokkos::sqrt(r[0] * r[0] + r[1] * r[1] + r[2] * r[2]);
    areaVec[0] = A[0];
    areaVec[1] = A[1];
    areaVec[2] = A[2];
    dv[0] = s * (r[0] - cx);
    dv[1] = s * (r[1] - cy);
    dv[2] = s * (r[2] - cz);
    conn[0] = r[0];
    conn[1] = r[1];
    conn[2] = r[2];
    return true;
  }
};
 
template <class Real, int MAXP, int MAXT, bool TrackAdj, class Policy = Voronoi>
KOKKOS_INLINE_FUNCTION void buildConvexCell(ConvexCell<Real, MAXP, MAXT, TrackAdj>& c,
                                            const Real L[3], const Real* relx, const Real* rely,
                                            const Real* relz, const int* ids, int nNbr,
                                            Real wSelf = Real(0), const Real* weights = nullptr) {
  c.initBox(L[0], L[1], L[2]);
  for (int i = 0; i < nNbr; ++i) {
    const Real n[3] = {relx[i], rely[i], relz[i]};
    Real wNbr = Real(0);
    if constexpr (Policy::kHasWeightDof) wNbr = weights[i];
    const Real off = Policy::template offsetFromRel<Real>(n, wSelf, wNbr);
    if constexpr (!Policy::kHasWeightDof) {
      if (!(off < Real(2) * c.maxVertexRsq()))
        break;  // security radius (bisector certificate; Voronoi only)
    }
    c.clip(n, off, ids[i]);
    if (c.overflow)
      return;
  }
}
 
}  // namespace peclet::voro
 
#endif  // PECLET_VORO_CONVEX_CELL_HPP