curvature.hpp

Go to the documentation of this file.

 
#ifndef PECLET_CORE_VOF_CURVATURE_HPP
#define PECLET_CORE_VOF_CURVATURE_HPP
 
#include <Kokkos_Core.hpp>
#include <Kokkos_MathematicalFunctions.hpp>
 
#include "peclet/core/vof/plic.hpp"
 
namespace peclet::core::vof {
 
// ---------------------------------------------------------------------------------------------
// cascade bookkeeping
// ---------------------------------------------------------------------------------------------
 
enum CurvatureBranch : int {
  kCurvNone = 0,        
  kCurvHf = 1,          
  kCurvHfMixed = 2,     
  kCurvHfFit = 3,       
  kCurvPv = 4,          
  kCurvPvReduced = 5,   
  kCurvNoEstimate = 6,  
};
 
inline constexpr int kHfColumn = 7;
inline constexpr int kPvHalf = 2;
inline constexpr double kPvWeightWidth = 2.5;
inline constexpr double kPtWeightWidth = 3.5;
 
inline constexpr double kHfPureEps = 1e-10;
 
KOKKOS_INLINE_FUNCTION bool hfIsFull(double c, double eps) {
  return c >= 1.0 - (eps > kHfPureEps ? eps : kHfPureEps);
}
KOKKOS_INLINE_FUNCTION bool hfIsEmpty(double c, double eps) {
  return c <= (eps > kHfPureEps ? eps : kHfPureEps);
}
 
// ---------------------------------------------------------------------------------------------
// tier 1 / 2 — height functions
// ---------------------------------------------------------------------------------------------
 
KOKKOS_INLINE_FUNCTION bool hfColumnHeight(const double* col, int nh, double& h, int& orient,
                                           double monoTol, double pureEps = kHfPureEps) {
  h = 0.0;
  orient = 0;
  const int R = nh / 2;
  double lo = 0.0, hi = 0.0;
  for (int k = 0; k < R; ++k) {
    lo += col[k];
    hi += col[nh - 1 - k];
  }
  if (lo > hi)
    orient = 1;
  else if (hi > lo)
    orient = -1;
  else
    return false;
 
  int a = -1, b = -1;
  if (orient == 1) {
    for (int k = R; k >= 0; --k)
      if (hfIsFull(col[k], pureEps)) {
        a = k;
        break;
      }
    for (int k = R; k < nh; ++k)
      if (hfIsEmpty(col[k], pureEps)) {
        b = k;
        break;
      }
  } else {
    for (int k = R; k >= 0; --k)
      if (hfIsEmpty(col[k], pureEps)) {
        a = k;
        break;
      }
    for (int k = R; k < nh; ++k)
      if (hfIsFull(col[k], pureEps)) {
        b = k;
        break;
      }
  }
  if (a < 0 || b < 0)
    return false;
 
  double s = 0.0;
  for (int k = a; k <= b; ++k) {
    s += col[k];
    // orient = +1 => non-increasing; orient = -1 => non-decreasing.
    if (k < b && (col[k + 1] - col[k]) * orient > monoTol)
      return false;
  }
  h = (orient == 1) ? (s + static_cast<double>(a - R) - 0.5)
                    : (s - static_cast<double>(b - R) - 0.5);
  return true;
}
 
KOKKOS_INLINE_FUNCTION double hfPatchKappa(const double h[9]) {
  const double hx = 0.5 * (h[2 + 3 * 1] - h[0 + 3 * 1]);
  const double hy = 0.5 * (h[1 + 3 * 2] - h[1 + 3 * 0]);
  const double hxx = h[2 + 3 * 1] - 2.0 * h[1 + 3 * 1] + h[0 + 3 * 1];
  const double hyy = h[1 + 3 * 2] - 2.0 * h[1 + 3 * 1] + h[1 + 3 * 0];
  const double hxy = 0.25 * (h[2 + 3 * 2] - h[0 + 3 * 2] - h[2 + 3 * 0] + h[0 + 3 * 0]);
  const double q = 1.0 + hx * hx + hy * hy;
  const double num = hxx + hyy + hxx * hy * hy + hyy * hx * hx - 2.0 * hxy * hx * hy;
  // kappa = 2H and H = -num / (2 q^{3/2}).
  return -num / (q * Kokkos::sqrt(q));
}
 
KOKKOS_INLINE_FUNCTION double hfPatchKappa(const double h[9], const VofMetric& g, int d) {
  const int d1 = (d + 1) % 3, d2 = (d + 2) % 3;
  const double r1 = g.h[d] / g.h[d1], r2 = g.h[d] / g.h[d2];
  const double hx = 0.5 * (h[2 + 3 * 1] - h[0 + 3 * 1]) * r1;
  const double hy = 0.5 * (h[1 + 3 * 2] - h[1 + 3 * 0]) * r2;
  const double hxx = (h[2 + 3 * 1] - 2.0 * h[1 + 3 * 1] + h[0 + 3 * 1]) * r1 / g.h[d1];
  const double hyy = (h[1 + 3 * 2] - 2.0 * h[1 + 3 * 1] + h[1 + 3 * 0]) * r2 / g.h[d2];
  const double hxy =
      0.25 * (h[2 + 3 * 2] - h[0 + 3 * 2] - h[2 + 3 * 0] + h[0 + 3 * 0]) * r1 / g.h[d2];
  const double q = 1.0 + hx * hx + hy * hy;
  const double num = hxx + hyy + hxx * hy * hy + hyy * hx * hx - 2.0 * hxy * hx * hy;
  return -num / (q * Kokkos::sqrt(q));
}
 
// ---------------------------------------------------------------------------------------------
// small linear algebra + geometry shared by the fit
// ---------------------------------------------------------------------------------------------
 
KOKKOS_INLINE_FUNCTION void curvFrame(const double n[3], double t1[3], double t2[3]) {
  int a = 0;
  double m = Kokkos::fabs(n[0]);
  if (Kokkos::fabs(n[1]) < m) {
    m = Kokkos::fabs(n[1]);
    a = 1;
  }
  if (Kokkos::fabs(n[2]) < m)
    a = 2;
  double e[3] = {0.0, 0.0, 0.0};
  e[a] = 1.0;
  // t1 = normalize(e x n)
  t1[0] = e[1] * n[2] - e[2] * n[1];
  t1[1] = e[2] * n[0] - e[0] * n[2];
  t1[2] = e[0] * n[1] - e[1] * n[0];
  const double s = Kokkos::sqrt(t1[0] * t1[0] + t1[1] * t1[1] + t1[2] * t1[2]);
  const double inv = (s > 0.0) ? 1.0 / s : 0.0;
  t1[0] *= inv;
  t1[1] *= inv;
  t1[2] *= inv;
  // t2 = n x t1
  t2[0] = n[1] * t1[2] - n[2] * t1[1];
  t2[1] = n[2] * t1[0] - n[0] * t1[2];
  t2[2] = n[0] * t1[1] - n[1] * t1[0];
}
 
KOKKOS_INLINE_FUNCTION double wendlandWeight(double r, double d) {
  if (!(d > 0.0))
    return 1.0;
  const double q = r / d;
  if (q >= 1.0)
    return 0.0;
  const double t = 1.0 - q;
  const double t2 = t * t;
  return (1.0 + 4.0 * q) * t2 * t2;
}
 
KOKKOS_INLINE_FUNCTION double paraboloidKappa(const double a[6]) {
  const double q = 1.0 + a[1] * a[1] + a[2] * a[2];
  const double num = a[5] * a[1] * a[1] - a[4] * a[1] * a[2] + a[3] * a[2] * a[2] + a[3] + a[5];
  return -2.0 * num / (q * Kokkos::sqrt(q));
}
 
KOKKOS_INLINE_FUNCTION bool curvSolveSym(const double Ain[6][6], const double bin[6], int n,
                                         double x[6], double pivTol) {
  double d[6];
  for (int i = 0; i < n; ++i) {
    if (!(Ain[i][i] > 0.0))
      return false;
    d[i] = 1.0 / Kokkos::sqrt(Ain[i][i]);
  }
  double A[6][6], b[6];
  for (int i = 0; i < n; ++i) {
    b[i] = bin[i] * d[i];
    for (int j = 0; j < n; ++j)
      A[i][j] = Ain[i][j] * d[i] * d[j];
  }
  // LDL^T in place (A[i][i] holds D, A[i][j] j<i holds L).
  for (int i = 0; i < n; ++i) {
    double s = A[i][i];
    for (int k = 0; k < i; ++k)
      s -= A[i][k] * A[i][k] * A[k][k];
    if (!(s > pivTol))  // NaN-safe: a NaN pivot fails the test
      return false;
    A[i][i] = s;
    for (int j = i + 1; j < n; ++j) {
      double t = A[j][i];
      for (int k = 0; k < i; ++k)
        t -= A[j][k] * A[i][k] * A[k][k];
      A[j][i] = t / s;
    }
  }
  // forward / diagonal / backward
  double y[6];
  for (int i = 0; i < n; ++i) {
    double s = b[i];
    for (int k = 0; k < i; ++k)
      s -= A[i][k] * y[k];
    y[i] = s;
  }
  for (int i = 0; i < n; ++i)
    y[i] /= A[i][i];
  for (int i = n - 1; i >= 0; --i) {
    double s = y[i];
    for (int k = i + 1; k < n; ++k)
      s -= A[k][i] * x[k];
    x[i] = s;
  }
  for (int i = 0; i < n; ++i)
    x[i] *= d[i];
  return true;
}
 
// ---------------------------------------------------------------------------------------------
// the PLIC polygon: plane ^ unit cube
// ---------------------------------------------------------------------------------------------
 
KOKKOS_INLINE_FUNCTION int plicPolygon(double mx, double my, double mz, double alpha,
                                       double v[8][3]) {
  const double m[3] = {mx, my, mz};
  const double nrm2 = m[0] * m[0] + m[1] * m[1] + m[2] * m[2];
  if (!(nrm2 > 0.0))
    return 0;
 
  int nv = 0;
  for (int q = 0; q < 8 && nv < 8; ++q) {
    const double A[3] = {static_cast<double>(q & 1), static_cast<double>((q >> 1) & 1),
                         static_cast<double>((q >> 2) & 1)};
    const double fA = m[0] * A[0] + m[1] * A[1] + m[2] * A[2] - alpha;
    for (int a = 0; a < 3 && nv < 8; ++a) {
      if ((q >> a) & 1)
        continue;  // walk each edge once, from its low corner
      double B[3] = {A[0], A[1], A[2]};
      B[a] = 1.0;
      const double fB = fA + m[a];
      // Straddle test with a half-open convention so a corner exactly on the plane is counted by
      // exactly one of its edges.
      const bool cross = (fA <= 0.0 && fB > 0.0) || (fB <= 0.0 && fA > 0.0);
      if (!cross)
        continue;
      const double den = fA - fB;
      const double t = (den != 0.0) ? fA / den : 0.0;
      double p[3] = {A[0], A[1], A[2]};
      p[a] = A[a] + t * (B[a] - A[a]);
      // de-duplicate (a corner hit is found by up to three edges)
      bool dup = false;
      for (int k = 0; k < nv; ++k) {
        const double dx = p[0] - v[k][0], dy = p[1] - v[k][1], dz = p[2] - v[k][2];
        if (dx * dx + dy * dy + dz * dz < 1e-24) {
          dup = true;
          break;
        }
      }
      if (dup)
        continue;
      v[nv][0] = p[0];
      v[nv][1] = p[1];
      v[nv][2] = p[2];
      ++nv;
    }
  }
  if (nv < 3)
    return nv;
 
  // Angular sort about the centroid in the plane's own frame.
  const double invn = 1.0 / Kokkos::sqrt(nrm2);
  const double nh[3] = {m[0] * invn, m[1] * invn, m[2] * invn};
  double e1[3], e2[3];
  curvFrame(nh, e1, e2);
  double cx = 0.0, cy = 0.0, cz = 0.0;
  for (int k = 0; k < nv; ++k) {
    cx += v[k][0];
    cy += v[k][1];
    cz += v[k][2];
  }
  const double inv = 1.0 / static_cast<double>(nv);
  cx *= inv;
  cy *= inv;
  cz *= inv;
  double ang[8];
  for (int k = 0; k < nv; ++k) {
    const double dx = v[k][0] - cx, dy = v[k][1] - cy, dz = v[k][2] - cz;
    ang[k] =
        Kokkos::atan2(dx * e2[0] + dy * e2[1] + dz * e2[2], dx * e1[0] + dy * e1[1] + dz * e1[2]);
  }
  for (int i = 1; i < nv; ++i) {  // insertion sort, nv <= 6
    const double a0 = ang[i];
    const double p0 = v[i][0], p1 = v[i][1], p2 = v[i][2];
    int j = i - 1;
    while (j >= 0 && ang[j] > a0) {
      ang[j + 1] = ang[j];
      v[j + 1][0] = v[j][0];
      v[j + 1][1] = v[j][1];
      v[j + 1][2] = v[j][2];
      --j;
    }
    ang[j + 1] = a0;
    v[j + 1][0] = p0;
    v[j + 1][1] = p1;
    v[j + 1][2] = p2;
  }
  return nv;
}
 
KOKKOS_INLINE_FUNCTION void polygonAreaCentroid(const double v[8][3], int nv, double ctr[3],
                                                double& area) {
  ctr[0] = ctr[1] = ctr[2] = 0.0;
  area = 0.0;
  if (nv < 3)
    return;
  for (int k = 1; k + 1 < nv; ++k) {
    const double a1[3] = {v[k][0] - v[0][0], v[k][1] - v[0][1], v[k][2] - v[0][2]};
    const double a2[3] = {v[k + 1][0] - v[0][0], v[k + 1][1] - v[0][1], v[k + 1][2] - v[0][2]};
    const double cx = a1[1] * a2[2] - a1[2] * a2[1];
    const double cy = a1[2] * a2[0] - a1[0] * a2[2];
    const double cz = a1[0] * a2[1] - a1[1] * a2[0];
    const double t = 0.5 * Kokkos::sqrt(cx * cx + cy * cy + cz * cz);
    area += t;
    ctr[0] += t * (v[0][0] + v[k][0] + v[k + 1][0]) / 3.0;
    ctr[1] += t * (v[0][1] + v[k][1] + v[k + 1][1]) / 3.0;
    ctr[2] += t * (v[0][2] + v[k][2] + v[k + 1][2]) / 3.0;
  }
  if (area > 0.0) {
    const double inv = 1.0 / area;
    ctr[0] *= inv;
    ctr[1] *= inv;
    ctr[2] *= inv;
  }
}
 
KOKKOS_INLINE_FUNCTION void polygonMoments2d(const double xy[8][2], int nv, double s[6]) {
  for (int i = 0; i < 6; ++i)
    s[i] = 0.0;
  if (nv < 3)
    return;
  for (int k = 0; k < nv; ++k) {
    const int k1 = (k + 1 == nv) ? 0 : k + 1;
    const double x0 = xy[k][0], y0 = xy[k][1];
    const double x1 = xy[k1][0], y1 = xy[k1][1];
    const double dx = x1 - x0, dy = y1 - y0;
    if (dy == 0.0)
      continue;
    s[0] += dy * 0.5 * (x0 + x1);
    s[1] += dy * (x0 * x0 + x0 * x1 + x1 * x1) / 6.0;
    s[2] += dy * (2.0 * x0 * y0 + x0 * y1 + x1 * y0 + 2.0 * x1 * y1) / 6.0;
    s[3] += dy * (x0 + x1) * (x0 * x0 + x1 * x1) / 12.0;
    // int_0^1 x^2 y dt, x = x0 + t dx, y = y0 + t dy
    s[4] += dy * 0.5 *
            (x0 * x0 * y0 + 0.5 * x0 * x0 * dy + x0 * dx * y0 + (2.0 / 3.0) * x0 * dx * dy +
             dx * dx * y0 / 3.0 + 0.25 * dx * dx * dy);
    // int_0^1 x y^2 dt
    s[5] += dy * (x0 * y0 * y0 + x0 * y0 * dy + x0 * dy * dy / 3.0 + 0.5 * dx * y0 * y0 +
                  (2.0 / 3.0) * dx * y0 * dy + 0.25 * dx * dy * dy);
  }
}
 
// ---------------------------------------------------------------------------------------------
// tier 2b — the mixed height function: a paraboloid through the interface positions that the
// consistent columns of ALL THREE directions do provide
// ---------------------------------------------------------------------------------------------
 
struct PtFit {
  double A[6][6];
  double b[6];
  int npt;
};
 
KOKKOS_INLINE_FUNCTION void ptFitInit(PtFit& f) {
  for (int i = 0; i < 6; ++i) {
    f.b[i] = 0.0;
    for (int j = 0; j < 6; ++j)
      f.A[i][j] = 0.0;
  }
  f.npt = 0;
}
 
KOKKOS_INLINE_FUNCTION void ptFitAdd(PtFit& f, const double X[3], const double org[3],
                                     const double t1[3], const double t2[3], const double nn[3],
                                     double dW, const VofMetric& g) {
  const double di[3] = {X[0] - org[0], X[1] - org[1], X[2] - org[2]};
  // Phase 3 (V2.4): the frame is orthonormal in PHYSICAL space, so the offsets are physical.
  // `g.h = {1,1,1}` makes this `di * 1.0`, i.e. today's numbers bit for bit.
  double d[3];
  g.toPhys(di, d);
  const double x = d[0] * t1[0] + d[1] * t1[1] + d[2] * t1[2];
  const double y = d[0] * t2[0] + d[1] * t2[1] + d[2] * t2[2];
  const double z = d[0] * nn[0] + d[1] * nn[1] + d[2] * nn[2];
  const double w = wendlandWeight(Kokkos::sqrt(x * x + y * y + z * z), dW);
  if (!(w > 0.0))
    return;
  const double phi[6] = {1.0, x, y, x * x, x * y, y * y};
  for (int i = 0; i < 6; ++i) {
    f.b[i] += w * phi[i] * z;
    for (int j = 0; j < 6; ++j)
      f.A[i][j] += w * phi[i] * phi[j];
  }
  ++f.npt;
}
 
KOKKOS_INLINE_FUNCTION void ptFitAdd(PtFit& f, const double X[3], const double org[3],
                                     const double t1[3], const double t2[3], const double nn[3],
                                     double dW) {
  ptFitAdd(f, X, org, t1, t2, nn, dW, VofMetric{});
}
 
KOKKOS_INLINE_FUNCTION bool ptFitSolve(const PtFit& f, double a[6], bool& red) {
  red = false;
  for (int i = 0; i < 6; ++i)
    a[i] = 0.0;
  if (f.npt >= 6 && curvSolveSym(f.A, f.b, 6, a, 1e-12))
    return true;
  if (f.npt < 3)
    return false;
  const int id[3] = {0, 3, 5};
  double Ar[6][6], br[6], ar[6];
  for (int i = 0; i < 3; ++i) {
    br[i] = f.b[id[i]];
    for (int j = 0; j < 3; ++j)
      Ar[i][j] = f.A[id[i]][id[j]];
  }
  if (!curvSolveSym(Ar, br, 3, ar, 1e-12))
    return false;
  a[0] = ar[0];
  a[3] = ar[1];
  a[5] = ar[2];
  red = true;
  return true;
}
 
// ---------------------------------------------------------------------------------------------
// tier 3 — the PLIC-volumetric paraboloid fit (Jibben 2019 / Han 2024 "PV")
// ---------------------------------------------------------------------------------------------
 
struct PvFit {
  double A[6][6];
  double b[6];
  int npoly;
};
 
KOKKOS_INLINE_FUNCTION void pvFitInit(PvFit& f) {
  for (int i = 0; i < 6; ++i) {
    f.b[i] = 0.0;
    for (int j = 0; j < 6; ++j)
      f.A[i][j] = 0.0;
  }
  f.npoly = 0;
}
 
struct PvTerm {
  double w, B, s[6];
  bool ok;
};
 
KOKKOS_INLINE_FUNCTION bool pvTermNormal(double mx, double my, double mz, const double t1[3],
                                         const double t2[3], const double nn[3], double cosMin,
                                         const VofMetric& g, double np[3]) {
  const double n2 = mx * mx + my * my + mz * mz;
  if (!(n2 > 0.0))
    return false;
  // Phase 3 (V2.4): the polygon's own normal enters the fit as a PHYSICAL direction (the frame is
  // physical), so it is pulled back through H^-1 and renormalized. At g.h = {1,1,1} the pullback
  // is `m/1.0` and `invn` is the same reciprocal square root as before — bitwise.
  const double mi[3] = {mx / g.h[0], my / g.h[1], mz / g.h[2]};
  const double invn = 1.0 / Kokkos::sqrt(mi[0] * mi[0] + mi[1] * mi[1] + mi[2] * mi[2]);
  np[0] = (mi[0] * t1[0] + mi[1] * t1[1] + mi[2] * t1[2]) * invn;
  np[1] = (mi[0] * t2[0] + mi[1] * t2[1] + mi[2] * t2[2]) * invn;
  np[2] = (mi[0] * nn[0] + mi[1] * nn[1] + mi[2] * nn[2]) * invn;
  if (!(np[2] > cosMin))
    return false;
  return true;
}
 
KOKKOS_INLINE_FUNCTION bool pvTermPolygon(PvTerm& t, const double np[3], const double v[8][3],
                                          int nv, const double off[3], const double org[3],
                                          const double t1[3], const double t2[3],
                                          const double nn[3], double dW, const VofMetric& g) {
  t.ok = false;
  // cell-local [0,1]^3 -> target-centred cell units -> the fit frame
  double xy[8][2];
  double zc = 0.0;
  double px = 0.0, py = 0.0;
  for (int k = 0; k < nv; ++k) {
    const double Xi[3] = {off[0] + v[k][0] - 0.5 - org[0], off[1] + v[k][1] - 0.5 - org[1],
                          off[2] + v[k][2] - 0.5 - org[2]};
    double X[3];
    g.toPhys(Xi, X);  // index displacement -> physical (identity at g.h = {1,1,1})
    xy[k][0] = X[0] * t1[0] + X[1] * t1[1] + X[2] * t1[2];
    xy[k][1] = X[0] * t2[0] + X[1] * t2[1] + X[2] * t2[2];
    const double z = X[0] * nn[0] + X[1] * nn[1] + X[2] * nn[2];
    px += xy[k][0];
    py += xy[k][1];
    zc += z;
  }
  const double invv = 1.0 / static_cast<double>(nv);
  px *= invv;
  py *= invv;
  zc *= invv;
 
  double s[6];
  polygonMoments2d(xy, nv, s);
  if (!(Kokkos::fabs(s[0]) > 1e-14))
    return false;
 
  // the polygon's own plane in the fit frame: z' = b0 + b1 x' + b2 y'
  const double b1 = -np[0] / np[2], b2 = -np[1] / np[2];
  const double b0 = zc - b1 * px - b2 * py;
 
  const double r = Kokkos::sqrt(px * px + py * py + zc * zc);
  const double w = wendlandWeight(r, dW);
  if (!(w > 0.0))
    return false;
 
  const double B = b0 * s[0] + b1 * s[1] + b2 * s[2];
  t.w = w;
  t.B = B;
  for (int i = 0; i < 6; ++i)
    t.s[i] = s[i];
  t.ok = true;
  return true;
}
 
KOKKOS_INLINE_FUNCTION bool pvFitTerm(PvTerm& t, double mx, double my, double mz, double alpha,
                                      const double off[3], const double org[3], const double t1[3],
                                      const double t2[3], const double nn[3], double dW,
                                      double cosMin, const VofMetric& g) {
  t.ok = false;
  double np[3];
  if (!pvTermNormal(mx, my, mz, t1, t2, nn, cosMin, g, np))
    return false;
  double v[8][3];
  const int nv = plicPolygon(mx, my, mz, alpha, v);
  if (nv < 3)
    return false;
  return pvTermPolygon(t, np, v, nv, off, org, t1, t2, nn, dW, g);
}
 
KOKKOS_INLINE_FUNCTION void pvFitAccum(PvFit& f, const PvTerm& t) {
  const double w = t.w, B = t.B;
  const double* s = t.s;
  for (int i = 0; i < 6; ++i) {
    f.b[i] += w * s[i] * B;
    for (int j = 0; j < 6; ++j)
      f.A[i][j] += w * s[i] * s[j];
  }
  ++f.npoly;
}
 
KOKKOS_INLINE_FUNCTION bool pvFitAdd(PvFit& f, double mx, double my, double mz, double alpha,
                                     const double off[3], const double org[3], const double t1[3],
                                     const double t2[3], const double nn[3], double dW,
                                     double cosMin, const VofMetric& g) {
  PvTerm t;
  if (!pvFitTerm(t, mx, my, mz, alpha, off, org, t1, t2, nn, dW, cosMin, g))
    return false;
  pvFitAccum(f, t);
  return true;
}
 
KOKKOS_INLINE_FUNCTION bool pvFitAdd(PvFit& f, double mx, double my, double mz, double alpha,
                                     const double off[3], const double org[3], const double t1[3],
                                     const double t2[3], const double nn[3], double dW,
                                     double cosMin) {
  return pvFitAdd(f, mx, my, mz, alpha, off, org, t1, t2, nn, dW, cosMin, VofMetric{});
}
 
KOKKOS_INLINE_FUNCTION bool pvFitSolve(const PvFit& f, double a[6], bool& red) {
  red = false;
  for (int i = 0; i < 6; ++i)
    a[i] = 0.0;
  if (f.npoly >= 6 && curvSolveSym(f.A, f.b, 6, a, 1e-12))
    return true;
  if (f.npoly < 3)
    return false;
  // reduced {1, x^2, y^2} == indices 0, 3, 5 of the full system
  const int id[3] = {0, 3, 5};
  double Ar[6][6], br[6], ar[6];
  for (int i = 0; i < 3; ++i) {
    br[i] = f.b[id[i]];
    for (int j = 0; j < 3; ++j)
      Ar[i][j] = f.A[id[i]][id[j]];
  }
  if (!curvSolveSym(Ar, br, 3, ar, 1e-12))
    return false;
  a[0] = ar[0];
  a[1] = 0.0;
  a[2] = 0.0;
  a[3] = ar[1];
  a[4] = 0.0;
  a[5] = ar[2];
  red = true;
  return true;
}
 
// ---------------------------------------------------------------------------------------------
// tier 3 at cost (design G, flow `doc/vof_curvature_cost_design.md` §5.1): per-cell caches of the
// PLIC polygon data, built ONCE per curvature pass, and the per-entry form of the accumulation
// ---------------------------------------------------------------------------------------------
 
KOKKOS_INLINE_FUNCTION void pvFitAccumLower(PvFit& f, const PvTerm& t) {
  const double w = t.w, B = t.B;
  const double* s = t.s;
  for (int i = 0; i < 6; ++i) {
    f.b[i] += w * s[i] * B;
    for (int j = 0; j <= i; ++j)
      f.A[i][j] += w * s[i] * s[j];
  }
  ++f.npoly;
}
 
inline constexpr int kPvEntries = 27;
 
KOKKOS_INLINE_FUNCTION void pvFitAccumEntry(double& a, const PvTerm& t, int e) {
  const double w = t.w;
  if (e >= 21) {
    const int i = e - 21;
    a += w * t.s[i] * t.B;
    return;
  }
  int i = 0;
  while ((i + 1) * (i + 2) / 2 <= e)
    ++i;
  const int j = e - i * (i + 1) / 2;
  a += w * t.s[i] * t.s[j];
}
 
KOKKOS_INLINE_FUNCTION void pvFitFromEntries(PvFit& f, const double ent[kPvEntries], int npoly) {
  for (int i = 0; i < 6; ++i) {
    f.b[i] = ent[21 + i];
    for (int j = 0; j < 6; ++j)
      f.A[i][j] = (j <= i) ? ent[i * (i + 1) / 2 + j] : 0.0;
  }
  f.npoly = npoly;
}
 
struct alignas(32) PvPolygon {
  double v[8][3];  
  double vbar[3];  
  int nv;          
};
 
KOKKOS_INLINE_FUNCTION void pvPolygonBuild(double mx, double my, double mz, double alpha,
                                           PvPolygon& P) {
  P.nv = plicPolygon(mx, my, mz, alpha, P.v);
  P.vbar[0] = P.vbar[1] = P.vbar[2] = 0.0;
  if (P.nv == 0)
    return;
  double sx = 0.0, sy = 0.0, sz = 0.0;
  for (int k = 0; k < P.nv; ++k) {
    sx += P.v[k][0];
    sy += P.v[k][1];
    sz += P.v[k][2];
  }
  const double inv = 1.0 / static_cast<double>(P.nv);
  P.vbar[0] = sx * inv - 0.5;
  P.vbar[1] = sy * inv - 0.5;
  P.vbar[2] = sz * inv - 0.5;
}
 
KOKKOS_INLINE_FUNCTION bool pvOutsideSupport(const PvPolygon& P, const double off[3],
                                             const double org[3], double dW, const VofMetric& g) {
  if (P.nv < 3)
    return false;
  if (!(dW > 0.0))
    return false;
  const double Xi[3] = {P.vbar[0] + off[0] - org[0], P.vbar[1] + off[1] - org[1],
                        P.vbar[2] + off[2] - org[2]};
  double X[3];
  g.toPhys(Xi, X);
  return X[0] * X[0] + X[1] * X[1] + X[2] * X[2] > dW * dW * (1.0 + 1e-9);
}
 
struct alignas(64) PvMoments {
  double n[3];   
  double c[3];   
  double a;      
  double m1[3];  
  double m2[6];  
};
 
KOKKOS_INLINE_FUNCTION void pvMomentsBuild(double mx, double my, double mz, double alpha,
                                           const VofMetric& g, PvMoments& P) {
  for (int a = 0; a < 3; ++a) {
    P.n[a] = 0.0;
    P.c[a] = 0.0;
    P.m1[a] = 0.0;
  }
  P.a = 0.0;
  for (int q = 0; q < 6; ++q)
    P.m2[q] = 0.0;
  const double m[3] = {mx, my, mz};
  if (vofPhysNormalInv(m, g, P.n) <= 0.0)
    return;
  double v[8][3];
  const int nv = plicPolygon(mx, my, mz, alpha, v);
  if (nv < 3)
    return;
  double Y[8][3];
  double sc[3] = {0.0, 0.0, 0.0};
  for (int k = 0; k < nv; ++k) {
    const double Xi[3] = {v[k][0] - 0.5, v[k][1] - 0.5, v[k][2] - 0.5};
    g.toPhys(Xi, Y[k]);
    sc[0] += Y[k][0];
    sc[1] += Y[k][1];
    sc[2] += Y[k][2];
  }
  const double inv = 1.0 / static_cast<double>(nv);
  P.c[0] = sc[0] * inv;
  P.c[1] = sc[1] * inv;
  P.c[2] = sc[2] * inv;
  constexpr int ip[6] = {0, 1, 2, 0, 0, 1}, iq[6] = {0, 1, 2, 1, 2, 2};
  for (int k = 1; k + 1 < nv; ++k) {
    const double* A = Y[0];
    const double* B = Y[k];
    const double* C = Y[k + 1];
    const double u[3] = {B[0] - A[0], B[1] - A[1], B[2] - A[2]};
    const double w[3] = {C[0] - A[0], C[1] - A[1], C[2] - A[2]};
    const double cx = u[1] * w[2] - u[2] * w[1];
    const double cy = u[2] * w[0] - u[0] * w[2];
    const double cz = u[0] * w[1] - u[1] * w[0];
    const double At = 0.5 * Kokkos::sqrt(cx * cx + cy * cy + cz * cz);
    const double S[3] = {A[0] + B[0] + C[0], A[1] + B[1] + C[1], A[2] + B[2] + C[2]};
    P.a += At;
    for (int a = 0; a < 3; ++a)
      P.m1[a] += At * S[a] / 3.0;
    for (int q = 0; q < 6; ++q) {
      const int p0 = ip[q], q0 = iq[q];
      P.m2[q] += At / 12.0 * (A[p0] * A[q0] + B[p0] * B[q0] + C[p0] * C[q0] + S[p0] * S[q0]);
    }
  }
}
 
KOKKOS_INLINE_FUNCTION bool pvFrameMoments(const PvMoments& P, const VofMetric& g, double nn[3],
                                           double t1[3], double t2[3], double org[3]) {
  if (!(P.n[0] * P.n[0] + P.n[1] * P.n[1] + P.n[2] * P.n[2] > 0.0))
    return false;
  nn[0] = P.n[0];
  nn[1] = P.n[1];
  nn[2] = P.n[2];
  curvFrame(nn, t1, t2);
  if (P.a > 0.0) {
    org[0] = P.m1[0] / P.a;
    org[1] = P.m1[1] / P.a;
    org[2] = P.m1[2] / P.a;
  } else {
    const double h[3] = {-0.5, -0.5, -0.5};
    g.toPhys(h, org);
  }
  return true;
}
 
KOKKOS_INLINE_FUNCTION bool pvTermMoments(PvTerm& t, const PvMoments& P, const double off[3],
                                          const double org[3], const double t1[3],
                                          const double t2[3], const double nn[3], double dW,
                                          double cosMin, const VofMetric& g) {
  t.ok = false;
  const double cj = P.n[0] * nn[0] + P.n[1] * nn[1] + P.n[2] * nn[2];
  if (!(cj > cosMin))  // also rejects a degenerate normal (n = 0) whenever cosMin >= 0
    return false;
  if (!(P.a > 0.0))
    return false;
  double d[3];
  g.toPhys(off, d);
  d[0] -= org[0];
  d[1] -= org[1];
  d[2] -= org[2];
  const double rc[3] = {P.c[0] + d[0], P.c[1] + d[1], P.c[2] + d[2]};
  const double r = Kokkos::sqrt(rc[0] * rc[0] + rc[1] * rc[1] + rc[2] * rc[2]);
  const double w = wendlandWeight(r, dW);
  if (!(w > 0.0))
    return false;
  const double a = P.a;
  const double* m1 = P.m1;
  const double* m2 = P.m2;
  const double Q1[3] = {m1[0] + d[0] * a, m1[1] + d[1] * a, m1[2] + d[2] * a};
  // Q2 = m2 + d m1^T + m1 d^T + d d^T a
  const double Qxx = m2[0] + 2.0 * d[0] * m1[0] + d[0] * d[0] * a;
  const double Qyy = m2[1] + 2.0 * d[1] * m1[1] + d[1] * d[1] * a;
  const double Qzz = m2[2] + 2.0 * d[2] * m1[2] + d[2] * d[2] * a;
  const double Qxy = m2[3] + d[0] * m1[1] + m1[0] * d[1] + d[0] * d[1] * a;
  const double Qxz = m2[4] + d[0] * m1[2] + m1[0] * d[2] + d[0] * d[2] * a;
  const double Qyz = m2[5] + d[1] * m1[2] + m1[1] * d[2] + d[1] * d[2] * a;
  const double u[3] = {Qxx * t1[0] + Qxy * t1[1] + Qxz * t1[2],
                       Qxy * t1[0] + Qyy * t1[1] + Qyz * t1[2],
                       Qxz * t1[0] + Qyz * t1[1] + Qzz * t1[2]};
  const double v[3] = {Qxx * t2[0] + Qxy * t2[1] + Qxz * t2[2],
                       Qxy * t2[0] + Qyy * t2[1] + Qyz * t2[2],
                       Qxz * t2[0] + Qyz * t2[1] + Qzz * t2[2]};
  double s[6];
  s[0] = cj * a;
  s[1] = cj * (t1[0] * Q1[0] + t1[1] * Q1[1] + t1[2] * Q1[2]);
  s[2] = cj * (t2[0] * Q1[0] + t2[1] * Q1[1] + t2[2] * Q1[2]);
  s[3] = cj * (t1[0] * u[0] + t1[1] * u[1] + t1[2] * u[2]);
  s[4] = cj * (t2[0] * u[0] + t2[1] * u[1] + t2[2] * u[2]);
  s[5] = cj * (t2[0] * v[0] + t2[1] * v[1] + t2[2] * v[2]);
  if (!(Kokkos::fabs(s[0]) > 1e-14))
    return false;
  t.w = w;
  t.B = cj * (nn[0] * Q1[0] + nn[1] * Q1[1] + nn[2] * Q1[2]);
  for (int i = 0; i < 6; ++i)
    t.s[i] = s[i];
  t.ok = true;
  return true;
}
 
}  // namespace peclet::core::vof
 
#endif  // PECLET_CORE_VOF_CURVATURE_HPP