core — VoF rung V3 (WO-O): interface curvature from volume fractions, as the Popinet (2009) height-function cascade with a PLIC-volumetric paraboloid-fit fallback. More...

#include <Kokkos_Core.hpp>
#include <Kokkos_MathematicalFunctions.hpp>
#include "peclet/core/vof/plic.hpp"

Classes
struct	peclet::core::vof::PtFit
	Accumulator for an ordinary weighted linear least-squares paraboloid through points. More...

struct	peclet::core::vof::PvFit
	Accumulator for the PV normal equations, Han et al. More...

struct	peclet::core::vof::PvTerm
	One stencil cell's PLIC polygon in the PV fit: the parameters of `pvFitTerm` and `pvFitAdd`. More...

struct	peclet::core::vof::PvPolygon
	The per-cell POLYGON cache entry (design G §5.1b, "V5"): `plicPolygon`'s output for one cell's PLIC plane, exactly as returned, plus the cell-centred INDEX-space vertex average used only by the support prefilter `pvOutsideSupport`. More...

struct	peclet::core::vof::PvMoments
	The per-cell MOMENT cache entry (design G §5.1d, "V6"): the 3-D area moments of a cell's PLIC polygon in PHYSICAL, cell-centred coordinates `Y = toPhys(v - 0.5)`. More...

Namespaces
namespace	peclet

namespace	peclet::core

namespace	peclet::core::vof

Enumerations
enum	peclet::core::vof::CurvatureBranch : int { peclet::core::vof::kCurvNone = 0 , peclet::core::vof::kCurvHf = 1 , peclet::core::vof::kCurvHfMixed = 2 , peclet::core::vof::kCurvHfFit = 3 , peclet::core::vof::kCurvPv = 4 , peclet::core::vof::kCurvPvReduced = 5 , peclet::core::vof::kCurvNoEstimate = 6 }
	Which branch of the cascade produced a cell's curvature. More...

Functions
bool	peclet::core::vof::hfIsFull (double c, double eps)
	Colour thresholds for "this cell is pure", at a caller-supplied tolerance.

bool	peclet::core::vof::hfIsEmpty (double c, double eps)

bool	peclet::core::vof::hfColumnHeight (const double *col, int nh, double &h, int &orient, double monoTol, double pureEps=kHfPureEps)
	One height from one column of colour values.

double	peclet::core::vof::hfPatchKappa (const double h[9])
	Mean-curvature-times-two from a 3x3 patch of heights, Han et al.

double	peclet::core::vof::hfPatchKappa (const double h[9], const VofMetric &g, int d)
	The same on an ANISOTROPIC cell (Phase 3; `flow/doc/anisotropic_vof.md` §4, decision V2.3).

void	peclet::core::vof::curvFrame (const double n[3], double t1[3], double t2[3])
	Right-handed orthonormal frame `(t1, t2, n)` from a unit vector `n`.

double	peclet::core::vof::wendlandWeight (double r, double d)
	Wendland C2 radial basis function, Han et al.

double	peclet::core::vof::paraboloidKappa (const double a[6])
	Mean curvature (times two) of the paraboloid `z = a0 + a1 x + a2 y + a3 x^2 + a4 x y + a5 y^2` at the origin — Han et al.

bool	peclet::core::vof::curvSolveSym (const double Ain[6][6], const double bin[6], int n, double x[6], double pivTol)
	Solve the symmetric positive-semidefinite `n x n` system `A x = b` (n <= 6) by LDL^T with Jacobi (diagonal) pre-scaling and a relative pivot guard.

int	peclet::core::vof::plicPolygon (double mx, double my, double mz, double alpha, double v[8][3])
	Vertices of the polygon `{x in [0,1]^3 : m .`

void	peclet::core::vof::polygonAreaCentroid (const double v[8][3], int nv, double ctr[3], double &area)
	Area and centroid of a planar polygon in 3-D (fan triangulation from `v[0]`).

void	peclet::core::vof::polygonMoments2d (const double xy[8][2], int nv, double s[6])
	The six monomial integrals `s = ( \|G\|, int x, int y, int x^2, int xy, int y^2 )` over the planar polygon `xy[0..nv)`, by Green's theorem — Han et al.

void	peclet::core::vof::ptFitInit (PtFit &f)

void	peclet::core::vof::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)
	Fold one interface point into the fit.

void	peclet::core::vof::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)
	Unit-metric overload (isotropic cells) — the pre-Phase-3 signature, unchanged arithmetic.

bool	peclet::core::vof::ptFitSolve (const PtFit &f, double a[6], bool &red)
	Solve the point fit, with the same rank cascade as `pvFitSolve`.

void	peclet::core::vof::pvFitInit (PvFit &f)

bool	peclet::core::vof::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])
	The NORMAL half of `pvFitTerm` (design G, flow `doc/vof_curvature_cost_design.md` §5.1a): the `n2 > 0` check, the pull-back `mi = m / g.h`, `invn = 1/\|mi\|`, the frame components `np = (mi.t1, mi.t2, mi.nn) * invn` and the `cosMin` guard — `pvFitTerm`'s expressions verbatim.

bool	peclet::core::vof::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)
	The POLYGON half of `pvFitTerm` (design G §5.1a): everything `pvFitTerm` does after `if (nv < 3) return false;`, verbatim — the vertex transform through `toPhys` into the fit frame, the vertex averages `px, py, zc`, `polygonMoments2d`, the `\|s0\| > 1e-14` test, the polygon's own plane `b0, b1, b2`, the Wendland radius and weight, and `B` — given the stencil cell's polygon `v[0..nv)` (`plicPolygon`'s output, cell-local [0,1]^3) and `np` from `pvTermNormal`.

bool	peclet::core::vof::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)
	One stencil polygon's contribution to the PV fit: everything `pvFitAdd` computes before the accumulation — the polygon, the frame transform, the moments `s`, the Wendland weight `w` and `B`.

void	peclet::core::vof::pvFitAccum (PvFit &f, const PvTerm &t)
	Fold one accepted term (`t.ok`) into the fit — the accumulation half of `pvFitAdd`, with its expressions verbatim.

bool	peclet::core::vof::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)
	Fold one stencil cell's PLIC polygon into the fit: `pvFitTerm`, then `pvFitAccum` if the polygon is accepted.

bool	peclet::core::vof::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)
	Unit-metric overload (isotropic cells) — the pre-Phase-3 signature, unchanged arithmetic.

bool	peclet::core::vof::pvFitSolve (const PvFit &f, double a[6], bool &red)
	Solve the accumulated PV system.

void	peclet::core::vof::pvFitAccumLower (PvFit &f, const PvTerm &t)
	Lower-triangle accumulation: `pvFitAccum` restricted to `j <= i` (design G §5.1c).

void	peclet::core::vof::pvFitAccumEntry (double &a, const PvTerm &t, int e)
	ONE entry `e` in [0, 27) of the accumulation, for an entry-parallel fold (design G §5.4: each entry summed serially in canonical order by its own lane): `e < 21` is `A[i][j]` with `e = i(i+1)/2 + j`, `j <= i`, and `a += t.w * t.s[i] * t.s[j]`; `e >= 21` is `b[e - 21]` and `a += t.w * t.s[e-21] * t.B`.

void	peclet::core::vof::pvFitFromEntries (PvFit &f, const double ent[kPvEntries], int npoly)
	A `PvFit` from the 27 entries of `pvFitAccumEntry` (+ `npoly`); the upper triangle is zero.

void	peclet::core::vof::pvPolygonBuild (double mx, double my, double mz, double alpha, PvPolygon &P)
	Build a `PvPolygon` from a cell's PLIC plane: `nv = plicPolygon(m, alpha, v)`, then `vbar`.

bool	peclet::core::vof::pvOutsideSupport (const PvPolygon &P, const double off[3], const double org[3], double dW, const VofMetric &g)
	true iff the cached polygon is CERTAINLY outside the Wendland support of a target with origin `org` (index units, as `pvFitTerm`'s), so `pvTermPolygon` would reject it with `w == 0`.

void	peclet::core::vof::pvMomentsBuild (double mx, double my, double mz, double alpha, const VofMetric &g, PvMoments &P)
	Build a `PvMoments` from a cell's PLIC plane (index-space `m`, `alpha`) and the cell metric: the physical unit normal, then — when the polygon has >= 3 vertices — its vertex average and its area moments by a fan from `Y_0`.

bool	peclet::core::vof::pvFrameMoments (const PvMoments &P, const VofMetric &g, double nn[3], double t1[3], double t2[3], double org[3])
	Tier 3's fit frame from the TARGET's cached moments: `nn = P.n`, `curvFrame(nn, t1, t2)`, and the origin `org` = the polygon's area centroid `m1 / a` in PHYSICAL target-centred units (`toPhys(-0.5, -0.5, -0.5)` when there is no polygon, today's degenerate `ctr = 0` case).

bool	peclet::core::vof::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)
	One stencil polygon's PV term from its cached moments (design G §5.1d, "V6"): the same integrals as `pvFitTerm`, transformed instead of re-integrated.

Variables
constexpr int	peclet::core::vof::kHfColumn = 7
	Column length of the height function, `NH` in Han et al.

constexpr int	peclet::core::vof::kPvHalf = 2
	Half-width of the PV fit stencil: S = 2*kPvHalf + 1 = 5, Han et al.'s recommended S.

constexpr double	peclet::core::vof::kPvWeightWidth = 2.5
	Wendland support width `d` in cell units.

constexpr double	peclet::core::vof::kPtWeightWidth = 3.5
	Wendland support width for the tier-2b height-position fit.

constexpr double	peclet::core::vof::kHfPureEps = 1e-10
	The FLOOR of the "this cell is pure" thresholds.

constexpr int	peclet::core::vof::kPvEntries = 27
	Number of independent entries of the PV normal equations: the 21 of the lower triangle (with the diagonal) and the 6 of `b`.

Detailed Description

core — VoF rung V3 (WO-O): interface curvature from volume fractions, as the Popinet (2009) height-function cascade with a PLIC-volumetric paraboloid-fit fallback.

Container-free by contract, exactly as vof/plic.hpp: every entry point takes scalars and small local arrays only — no Kokkos::View, no grid indexing, no halo/topology types. These kernels are scheduled for promotion to peclet::core::vof at the V4 freeze (suite/docs/VOF_PLAN.md §11), where the structured solver, the AMR octree and the per-bubble block container share ONE copy, and the promotion must be a plain file move. The view-level driver that walks the block, owns the scratch fields and registers "kappa" lives next door in vof/curvature_field.hpp — the same split plic.hpp / colour_field.hpp already uses.

The cascade

Standard height function in the preferred direction (the Cartesian axis with the largest |n_d| of the cell's own PLIC normal). Heights are column sums of C over NH = 7 cells — Han/Evrard/Desjardins (2024) eq. (3), the column length used throughout the literature — over the 3x3 transverse patch, then the second-order finite differences of their eqs. (4)-(5).
Mixed-direction HF, in two forms, of which only the first ships on:
- 2a the same construction in the two remaining directions, tried in order of decreasing |n_d|. This is what height_curvature() in Basilisk's curvature.h does (foreach_dimension() over the sorted normal), and it is free — the columns are the same colour data read along another axis. Measured: it fires on 0.00-0.06 % of a sphere's interfacial cells, because the direction with the largest |n_d| is by construction the one most likely to close.
- 2b Popinet's generalised height function / Basilisk's height_curvature_fit: a paraboloid through the interface positions of whichever of the 27 columns DO close. Implemented (PtFit) and default OFF — it is measurably worse than routing those cells to tier 3. See VofCurvature::useMixedHeightFit for the numbers and the mechanism, and the DEVIATION note at the bottom of this comment for its halo-imposed restriction.
PLIC-volumetric paraboloid fit (the "PV" method of Jibben et al. 2015, as evaluated and recommended by Han, Evrard & Desjardins, IJMF 174:104769 (2024), arXiv:2304.08643) on a 5^3 stencil with a Wendland C2 radial weight of width d = 2.5 cells — their measured best cost/accuracy fallback, and the exact configuration their §5 dynamic tests use.

Why curvature never touches a normal derivative. Rung V0 measured that MYC's normal error does not converge under refinement (order 0.83; ~28 % of mixed cells take the non-converging Youngs branch) while its reconstruction error does (1.98) — doc/vof_workorders.md, WO-D finding 2. Every branch here therefore takes its geometry from column sums of C (tiers 1-2) or from the PLIC planes' volumes (tier 3). The MYC normal appears in exactly two roles, both categorical rather than differentiated: it orders the three candidate column directions (Basilisk does the same), and it defines the local frame of the paraboloid fit. Neither is differenced, and an O(1 deg) error in either changes only which equally-valid representation of the same interface is fitted.

Sign and units

Curvature is returned as kappa = 2H in units of 1/h (cell units; multiply by 1/h for physical units), positive for a convex blob of liquid — a sphere of liquid of radius R cells gives kappa = +2/R. C is the liquid fraction and the PLIC normal m points into the gas (plic.hpp convention), so the two conventions are consistent by construction.

The height is a signed column sum (hfColumnHeight), and the sign is what makes one formula serve both orientations: with liquid on the low side the sum equals z_I + R + 1/2, with liquid on the high side it equals R + 1/2 - z_I, so in both cases h grows with the amount of liquid — and Han's eq. (5), whose numerator is odd in h and whose denominator is even, returns the same H either way. Only the consistency test cares about the orientation: all nine columns of a patch must be oriented the same way, or the nine numbers do not describe one single-valued surface. (h itself is the interface position along the column only up to that sign: the position is orient * h, which is what tier 2b's point cloud needs.)

Two facts from the literature the callers must respect, not fight

A cascade that never falls back has a bug. Han et al. measure the HF method failing for ~51 % of random paraboloids with curvedness C > 1e-1, and — on the 2D stationary droplet, i.e. the friendliest possible case — still up to 0.9 % of interfacial cells at D/dx = 102.4. Below ~4-5 cells per diameter it fails everywhere.
Curvature error stops converging below C*dx ~ 1e-2 with advection-realistic volume fractions, for every published method (their §3 and §4.1, following Remmerswaal & Veldman's Lemma 3: a first-order-in-Linf volume fraction field makes the curvature zeroth-order). A plateau on fine grids is the method's physics; it is measured and reported, never chased.

DEVIATION recorded: which "mixed-direction HF" this is, and why 2b ships off

Popinet (2009) §3 and Basilisk's curvature.h offer two things that both get called "mixed": (a) trying the three column directions in turn (height_curvature), and (b) a paraboloid fit to the interface positions implied by every defined height, in any direction (height_curvature_fit). Both are implemented here as tiers 2a and 2b.

2b carries one restriction and one measured verdict.

Restriction. Basilisk gathers heights over the 3^3 neighbourhood, which needs colour at +/-1 +/- 3 = +/-4 cells, i.e. a g = 4 halo or a separate exchange of the three height fields. This implementation gathers from the 3x3 patches of the three directions instead — cells at offsets (+/-1, +/-1, +/-3) — so the WHOLE rung fits exactly inside the colour field's existing g = 3 halo (tier 1 reaches 3x3x7, tier 3 reaches a 5^3 of MYC stencils = 7^3 of colour; both are exactly +/-3). No new halo machinery is needed and none was added.
Verdict: 2b ON destroys the convergence of the MAX curvature error (order 0.00 vs 1.86 over 16^3 -> 64^3 on the exact-fraction sphere), on the very cells tier 3 handles at second order, for every Wendland width from 1.5 to 6.0 cells. The mechanism is structural — its data set is the columns the height function could close, a slope-selected and therefore asymmetric subset, whose lever-arm bias does not shrink with h. Full numbers and the reasoning in VofCurvature::useMixedHeightFit; regenerated on every ctest run by gate G of tests/kokkos/test_vof_curvature.cpp. Han et al. independently measure the PV fallback as better than HF once volume fractions carry transport error — which is exactly the regime the fallback fires in.

Sources

S. Popinet, An accurate adaptive solver for surface-tension-driven interfacial flows, JCP 228:5838 (2009) — the HF cascade and the generalised height function.
A. Han, F. Evrard, O. Desjardins, Comparison of methods for curvature estimation from volume fractions, IJMF 174:104769 (2024), arXiv:2304.08643 — eqs. (3)-(5) (HF), (6), (8) (Wendland), (10) (curvature of a paraboloid), (11)-(15) (the PV normal equations), and every quantitative claim quoted above.
Z. Jibben, N. N. Carlson, M. M. Francois, A paraboloid fitting technique for calculating curvature from piecewise-linear interface reconstructions on 3D unstructured meshes, Comput. Math. Appl. (2019), arXiv:1712.05467 — the PV method itself.
H. Wendland, Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree, Adv. Comput. Math. 4:389 (1995) — the C2 kernel of eq. (8).

**Layer L1 (VOF_PLAN.md §11): promoted from flow/src/vof/ to core (WO-W0, 2026-09-02).** These are the container-free VoF math kernels — one copy shared by every container: flow's structured colour field, the per-bubble block container (Part III) and the future AMR path. flow/src/vof/curvature.hpp is now a thin include + a using-directive, so every existing peclet::flow::vof:: spelling keeps working unchanged.

Definition in file curvature.hpp.

Classes

Namespaces

Enumerations