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

struct	PvFit
	Accumulator for the PV normal equations, Han et al. More...

struct	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...

struct	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	PvTerm
	One stencil cell's PLIC polygon in the PV fit: the parameters of `pvFitTerm` and `pvFitAdd`. More...

struct	VofMetric
	The per-axis cell size, in the solver's reference length (`hRef`), of the cell these kernels map to the unit cube. More...

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

enum	VofCellKind : unsigned char { kVofFluid = 0 , kVofSolid = 1 }
	Cell classification on the colour block. More...

enum	VofFillState : unsigned char { kVofFillNone = 0 , kVofFillFluid = 1 , kVofFillPass0 = 1 }
	Band-fill state, stored per cell on the colour block. More...

enum	VofWettingPivot : int { kVofPivotVolume = 0 , kVofPivotInterface = 1 , kVofPivotWallNormal = 2 , kVofPivotContactLine = 3 }
	How the theta-plane is anchored in the fluid cell. More...

enum	VofWettingBranch : int { kVofWetNone = 0 , kVofWetTheta = 1 , kVofWetNeighbour = 2 , kVofWetPure = 3 , kVofWetParallel = 4 , kVofWetNeutral = 5 , kVofWetCount = 6 }
	Which rule produced a band cell's colour, for `contact_angle_diagnostics()`. More...

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

bool	hfIsEmpty (double c, double eps)

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

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

double	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	curvFrame (const double n[3], double t1[3], double t2[3])
	Right-handed orthonormal frame `(t1, t2, n)` from a unit vector `n`.

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

double	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	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	plicPolygon (double mx, double my, double mz, double alpha, double v[8][3])
	Vertices of the polygon `{x in [0,1]^3 : m .`

void	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	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	ptFitInit (PtFit &f)

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)
	Fold one interface point into the fit.

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

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

void	pvFitInit (PvFit &f)

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])
	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	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	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	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	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	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	pvFitSolve (const PvFit &f, double a[6], bool &red)
	Solve the accumulated PV system.

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

void	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	pvFitFromEntries (PvFit &f, const double ent[kPvEntries], int npoly)
	A `PvFit` from the 27 entries of `pvFitAccumEntry` (+ `npoly`); the upper triangle is zero.

void	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	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	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	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	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.

double	vofEpsEff (double eps, bool solid)
	The effective fluid volume of a cell (rule 1).

bool	vofIsSolidCell (double eps, const double oLo[3], const double oHi[3])
	Is this cell SOLID? `eps == 0` and every one of the six faces closed.

double	vofCutFlux (double openness, double rawFlux)
	The openness-weighted geometric flux (1): one number per face, used with opposite signs by the two neighbours, so the sum telescopes bit-exactly whatever `o` is.

double	vofCutFluxClamp (double fluxMag, double sweep, double epsDon, double cDon)
	Weymouth's admissible flux interval, generalized to a cut donor cell — the cut-cell twin of `momentum_advect.hpp` point 3, and for the same reason.

double	vofCutCourant (double openness, double a, double eps)
	The boundedness-relevant Courant number of a face against a cell (rule 2).

double	vofClipCut (double c, double epsEff, double &moved)
	The clip of rule 3, returned as the CLIPPED value; `moved` receives the liquid VOLUME the clip created (positive) or destroyed (negative) in this cell, i.e.

bool	vofFillReadable (unsigned char state, int k)
	Is `state` readable by pass `k` (1-based)? Fluid always; a solid cell only if an EARLIER pass filled it.

int	vofGhostDepth1 (int i, int g, int n)
	Ghost depth of an index on one axis of an extended block: 0 for an inner cell, growing outward.

int	vofGhostDepth (int x, int y, int z, int g, int nx, int ny, int nz)
	Ghost depth of a cell = the max over the three axes.

double	plicSq (double x)

double	plicCube (double x)

double	plicNormalizeL1 (double m[3])
	L1-normalize a normal in place. Returns the L1 norm before normalization (0 => untouched).

double	vofPhysNormal (const double m[3], const VofMetric &g, double n[3])
	The PHYSICAL unit normal of an index-space normal `m` (any scale), i.e.

double	vofPhysNormalInv (const double m[3], const VofMetric &g, double n[3])
	`vofPhysNormal` in the RECIPROCAL-MULTIPLY form — `n = (H^-1 m) * (1/\|H^-1 m\|)`, which is what the curvature cascade's and the interface area's frame construction has always computed.

void	vofIndexNormal (const double n[3], const VofMetric &g, double m[3])
	The INDEX-space normal of a physical normal `n` (any scale): `m = H n`, returned unnormalized (`plicAlpha` / `plicVolume` renormalize internally, so no L1 pass is needed).

double	plicVolumeReduced (double w, double n1, double n2, double n3)
	Reduced-symmetry forward relation.

double	plicVolume (double mx, double my, double mz, double alpha)
	Volume fraction of the unit cell [0,1]^3 lying on the fluid side `m .`

double	plicAlphaReduced (double v, double n1, double n2, double n3)
	Reduced-symmetry analytic inverse of `plicVolumeReduced`.

double	plicAlpha (double mx, double my, double mz, double v)
	Plane offset alpha such that `plicVolume(mx,my,mz,alpha) == V`, exactly (analytic inverse).

double	plicSlabVolume (double mx, double my, double mz, double alpha, int dir, double a, double b)
	Fluid volume of the PLIC polyhedron inside the axis-aligned slab `a <= x_dir <= b`, expressed as a fraction of the WHOLE cell (so `plicSlabVolume(..., 0, 1) == plicVolume(...)`).

double	plicBoxVolume (double mx, double my, double mz, double alpha, double a0, double b0, double a1, double b1, double a2, double b2)
	Fluid volume of the PLIC polyhedron inside the axis-aligned BOX `[a0,b0] x [a1,b1] x [a2,b2]`, as a fraction of the WHOLE cell.

double	faceFluxVolume (double mx, double my, double mz, double alpha, int dir, double f)
	Donor-cell geometric flux volume: the fluid volume swept out of the slab `[0, f]` along `dir`, as a fraction of the whole cell.

int	plicSt (int i, int j, int k)
	Stencil index: x-fastest, offsets (i-1, j-1, k-1) with i,j,k in {0,1,2} (`suite/docs/CONVENTIONS.md`).

void	youngsNormal (const double c[27], double m[3])
	Youngs' normal: the 27-point weighted gradient (corner 1, edge 2, face 4 on each of the two opposing planes), L1-normalized, pointing into the gas (m ~ -grad C).

void	mycNormal (const double c[27], double m[3])
	MYC — Mixed Youngs-Centred normal (Aulisa et al.

double	planeCellFraction (double mx, double my, double mz, double alphaGlobal, double x0, double y0, double z0, double h)
	Exact fluid fraction of the axis-aligned cell `[x0,x0+h] x [y0,y0+h] x [z0,z0+h]` under the GLOBAL plane `m .`

double	planeCellFractionAniso (double mx, double my, double mz, double alphaGlobal, double x0, double y0, double z0, const VofMetric &g)
	The same for an anisotropic cell `[x0, x0+g.h[0]] x [y0, y0+g.h[1]] x [z0, z0+g.h[2]]`: the index normal of the global plane is `H n` (`vofIndexNormal`), so the exact fraction is one `plicVolume` call.

double	sphereCellFraction (double cx, double cy, double cz, double r, double x0, double y0, double z0, double h, int levels)
	Fluid fraction of a cell for a sphere (fluid = INSIDE the sphere), by recursive octree subdivision with an exact tangent-plane closure at the leaves.

bool	youngsNormalFluidOnly (const double c[27], const unsigned char fluid[27], double m[3])
	Youngs' 27-point weighted gradient restricted to non-solid cells.

void	vofPlicCentroid (const double m[3], double alpha, double p[3])
	The centroid of the PLIC polygon of cell `f` in the cell's own unit-cube frame, for the plane `m .`

int	vofWettingPlane (const double mfIn[3], double cf, const double nwIn[3], double cosT, double sinT, double sdfF, int pivot, double tEps, double mth[3], double &alphaTh, double &cosApp, const VofMetric &g)
	Build the theta-plane in the FLUID cell's own unit-cube frame.

int	vofWettingPlane (const double mfIn[3], double cf, const double nwIn[3], double cosT, double sinT, double sdfF, int pivot, double tEps, double mth[3], double &alphaTh, double &cosApp)
	Unit-metric overload (isotropic cells) — the pre-Phase-3 signature, unchanged arithmetic.

double	vofWettingFraction (const double mth[3], double alphaTh, const int ds[3])
	The liquid fraction of the solid cell at integer offset `ds` from the fluid cell, under the theta-plane returned by `vofWettingPlane`.

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

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

constexpr double	kPvWeightWidth = 2.5
	Wendland support width `d` in cell units.

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

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

constexpr int	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`.

constexpr double	kVofEpsFloor = 1.0 / 64.0
	The smallest fluid volume `buildCellFraction`'s 4³ subsampling can resolve (one subcell of 64).

constexpr double	kVofCourantEpsFloor = 0.1
	Below this the interface-local Courant number stops being amplified by 1/eps (rule 2): a cut cell at eps = 0.1 already reports 10x its face Courant number, and chasing thinner slivers with dt buys boundedness the clip provides for free.

Enumeration Type Documentation

◆ CurvatureBranch

enum peclet::core::vof::CurvatureBranch : int

Which branch of the cascade produced a cell's curvature.

Written to a companion field so the branch census is a measurement and never an inference, and so a cell that got NO estimate is LOUD (kCurvNoEstimate) instead of carrying a silently-zero curvature.

Enumerator
kCurvNone	not an interfacial cell — kappa is 0 because there is no interface
kCurvHf	tier 1: standard HF in the preferred (largest \|n_d\|) direction
kCurvHfMixed	tier 2a: standard HF in one of the two remaining directions
kCurvHfFit	tier 2b: paraboloid fit to the interface positions of whichever columns ARE consistent, pooled over all three directions
kCurvPv	tier 3: PLIC-volumetric paraboloid fit, 5^3 Wendland-weighted
kCurvPvReduced	tier 3 with the rank-deficient 3-parameter model (a1 = a2 = a4 = 0)
kCurvNoEstimate	NO estimate: fewer than 3 usable PLIC polygons in the stencil.

Definition at line 134 of file curvature.hpp.

◆ VofCellKind

enum peclet::core::vof::VofCellKind : unsigned char

Cell classification on the colour block.

Kept as an unsigned char field so it can ride the ordinary halo exchange as a double and be read as a flag in every kernel.

Enumerator
kVofFluid	carries colour; `eps_eff > 0`; at least one open face OR `eps > 0`
kVofSolid	`eps == 0` AND all six faces closed: no colour, filled for the stencils

Definition at line 94 of file cutcell.hpp.

◆ VofFillState

enum peclet::core::vof::VofFillState : unsigned char

Band-fill state, stored per cell on the colour block.

kVofFillFluid marks a cell whose colour is data (a fluid cell); kVofFillPass0 + k marks a solid cell filled in pass k (k = 1..3); kVofFillNone is an unfilled solid cell. Pass k reads a neighbour iff its state is kVofFillFluid or kVofFillPass0 + j with j < k — never a cell being written in the same pass, which is what makes the three passes race-free without a second buffer.

Enumerator
kVofFillNone
kVofFillFluid
kVofFillPass0	pass k marks `kVofFillPass0 + k`, i.e. 2, 3, 4

Definition at line 195 of file cutcell.hpp.

◆ VofWettingPivot

enum peclet::core::vof::VofWettingPivot : int

How the theta-plane is anchored in the fluid cell.

See the file header for why the default is not the work order's rule.

Enumerator
kVofPivotVolume	`plicAlpha(m_theta, C_f)`: volume-consistent. DEFAULT, idempotent.
kVofPivotInterface	through `p_f` (Afkhami-Bussmann / Basilisk). Idempotent.
kVofPivotWallNormal	WO-S as written: `c = p_f - sdf(p_f) n_w`. NOT idempotent.
kVofPivotContactLine	the point of the fluid plane on the wall nearest `p_f`. Idempotent.

Definition at line 126 of file wetting.hpp.

◆ VofWettingBranch

enum peclet::core::vof::VofWettingBranch : int

Which rule produced a band cell's colour, for contact_angle_diagnostics().

Enumerator
kVofWetNone	not written by pass 1
kVofWetTheta	the theta-plane of this column's own anchor cell
kVofWetNeighbour	the mean of the theta-planes of the anchor's MIXED neighbours, used where the anchor itself is pure phase but the contact line passes through the column next to it — see "the anchor's own column is not enough" in the file header
kVofWetPure	pure-phase continuation (`C_f` is 0 or 1, no interface in reach)
kVofWetParallel	interface parallel to the wall (\|t\| < tEps): `m_f` kept unrotated
kVofWetNeutral	fell back to WO-Q's neutral mean (no usable `n_w`/`f`/`m_f`)
kVofWetCount

Definition at line 134 of file wetting.hpp.

Function Documentation

◆ hfIsFull()

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

Colour thresholds for "this cell is pure", at a caller-supplied tolerance.

Why this is a parameter and not the constant it used to be. The height function is the only consumer in the VoF stack that decides what a PURE cell is on its own. Every other one is told: the advector's wispEps, phase change's pcEffInterfaceEps(). When a caller runs the advector at a looser tolerance than this floor – enable_vof sets wispEps = 1e-8 – the colour field legitimately carries bulk cells at 1 - O(1e-9): the advector calls them pure, so it fluxes them algebraically and never reconstructs them back onto exactly 1.0. Judging those same cells at 1e-10 makes the column walk fail to find a pure end, and the whole height-function tier collapses into the paraboloid fallback – silently, because a fallback is not an error. That is measured: on the mode-2 droplet the HF tier fell from 790 cells to 134 (37 % -> 89 % fallback) over 2.5 periods, tracking 425 bulk cells drifting into the band between the two tolerances.

The residue argument is unchanged and covers the whole range a caller may ask for: a cell accepted as full at 1 - eps contributes eps to the column sum, i.e. eps cells of height error, which at the 1e-8 the advector uses is eight orders below the height's own O(h^2) error.

eps <= kHfPureEps reproduces the pre-parameter arithmetic bit for bit, which is what the default argument on hfColumnHeight preserves for every caller that does not pass one.

Definition at line 185 of file curvature.hpp.

References kHfPureEps.

Referenced by hfColumnHeight().

◆ hfIsEmpty()

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

Definition at line 188 of file curvature.hpp.

References kHfPureEps.

Referenced by hfColumnHeight().

◆ hfColumnHeight()

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.

Parameters

col	`nh` colour values along the column, `col[nh/2]` being the target cell.
nh	column length (odd; `kHfColumn`).
h	[out] the height in cell units, relative to the centre cell's centre. Signed so that it increases with liquid whichever end the liquid is at.
orient	[out] +1 when the low end of the column is liquid, -1 when the high end is.
monoTol	tolerance on the monotonicity test (see below).
pureEps	tolerance of the PURITY test that finds the window ends (`hfIsFull`/`hfIsEmpty`). Pass the tolerance the COLOUR FIELD was advected at – `WyAdvector::wispEps` – so that the cells the advector treats as pure are the cells this walk accepts as a window end. Floored at `kHfPureEps`; the default reproduces the pre-parameter arithmetic exactly.

Returns: true iff the column is CONSISTENT.

This is Popinet's outward accumulation, not a fixed symmetric window: starting at the centre cell, walk down to the first pure cell and up to the first pure cell of the opposite kind, and integrate over that window [a, b] only. The height then follows from the window's own end:

liquid low  (col[a] full,  col[b] empty):  h =  sum_{a..b} C + a - 1/2
liquid high (col[a] empty, col[b] full ):  h =  sum_{a..b} C - b - 1/2

(indices centre-relative). When both ends of the full column happen to be pure the two forms collapse to the familiar sum(col) - NH/2 of Han et al. eq. (3) — the window formulation is strictly more permissive, because it also serves a column whose far end is not pure (a second interface just outside, a wisp) and a column whose centre cell is itself pure (which every neighbour column of a sloping patch eventually is).

A column is consistent when such a window exists AND the colour is monotone across it to within monoTol. Monotonicity is what the work order's "where monotone columns exist" names, and it is not decoration: [1,1,0,1,0,0,0] brackets a transition but its sum is not a height, it is a height plus a sub-grid blob.

The orientation is decided by which half of the column carries more liquid, so it is a property of the data rather than of an externally supplied normal — the cascade never differentiates a normal (see the file header). An exactly balanced column (a symmetric double interface) is rejected.

Definition at line 233 of file curvature.hpp.

References hfIsEmpty(), and hfIsFull().

◆ hfPatchKappa() [1/2]

double peclet::core::vof::hfPatchKappa ( const double h[9] )

Mean-curvature-times-two from a 3x3 patch of heights, Han et al.

eqs. (4)-(5) with dx = dy = 1 (cell units).

Parameters

h	nine heights, indexed `p + 3*q` with p, q in {0,1,2} standing for the transverse offsets {-1, 0, +1}. Which of the two transverse axes is p and which is q does not matter: eq. (5) is symmetric under the swap.

Returns: kappa = 2H, in 1/cell.

Definition at line 296 of file curvature.hpp.

◆ hfPatchKappa() [2/2]

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).

hfColumnHeight returns the height in CELLS along the column direction d, so the graph of the interface over the transverse plane is the PHYSICAL function f = h_d * h_idx of the physical transverse coordinates (h_d1 x_idx, h_d2 y_idx). Han et al. eqs. (4)-(5) on those spacings are the isotropic expressions times five per-axis ratios:

f_x  = hx * (h_d/h_d1)          f_xx = hxx * (h_d/h_d1^2)
f_y  = hy * (h_d/h_d2)          f_yy = hyy * (h_d/h_d2^2)
                                f_xy = hxy * (h_d/(h_d1 h_d2))

and the returned curvature is in 1/hRef — the unit ‘sigma’ kappa' = p'` already assumes.

Each ratio is a quotient of two entries of g.h, so at equal spacings every one of them is exactly 1.0 and hx * 1.0 == hx bitwise: the isotropic call reproduces hfPatchKappa(h) to the last bit (the one-argument overload above is what the unit-metric path still calls).

Parameters

d	the column direction (0/1/2); the transverse axes are `d1 = (d+1)%3`, `d2 = (d+2)%3`, in the order `curvHeightCell` fills the patch.

Definition at line 327 of file curvature.hpp.

References peclet::core::vof::VofMetric::h.

◆ curvFrame()

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.

Deterministic (the seed axis is the one with the smallest |n_a|), so the frame — and therefore every fitted coefficient — is a pure function of n and reproduces bitwise on any backend and any decomposition.

Definition at line 348 of file curvature.hpp.

Referenced by plicPolygon(), and pvFrameMoments().

◆ wendlandWeight()

double peclet::core::vof::wendlandWeight	(	double	r,
		double	d
	)

Wendland C2 radial basis function, Han et al.

eq. (8): (1 + 4 r/d)(1 - r/d)^4 on 0 <= r <= d, zero beyond.

Definition at line 376 of file curvature.hpp.

Referenced by ptFitAdd(), pvTermMoments(), and pvTermPolygon().

◆ paraboloidKappa()

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.

eq. (10), doubled. Sign convention: the paraboloid's z axis is the interface normal pointing into the GAS, so a liquid sphere of radius R gives +2/R.

Definition at line 390 of file curvature.hpp.

◆ curvSolveSym()

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.

Returns false — leaving x untouched — when the system is rank deficient, which is the caller's cue to drop to a smaller model rather than to emit a garbage curvature.

A is taken by value in the caller's local storage; this routine overwrites its own copy only.

Definition at line 402 of file curvature.hpp.

Referenced by ptFitSolve(), and pvFitSolve().

◆ plicPolygon()

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 .

x = alpha}, ordered around the polygon.

A plane cuts a cube in at most a hexagon, so v needs 6 slots (8 given, for the degenerate corner-coincidence cases before de-duplication). Edge walk + angular sort in the plane's own 2-D frame — the standard construction, and the only piece of polygon machinery this rung needs (the PV normal equations are line integrals around exactly this polygon, Han et al. eq. (14)). plic.hpp deliberately has no clipper; this one is confined to the curvature fallback.

Returns: the vertex count (0 when the plane misses the cube or m is degenerate).

Definition at line 465 of file curvature.hpp.

References peclet::core::cross(), and curvFrame().

Referenced by pvFitTerm(), pvMomentsBuild(), pvPolygonBuild(), and vofPlicCentroid().

◆ polygonAreaCentroid()

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]).

Definition at line 553 of file curvature.hpp.

Referenced by vofPlicCentroid().

◆ polygonMoments2d()

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.

eqs. (14a)-(14f).

DEVIATION recorded. These are re-derived here rather than transcribed, using int_G x^a y^b dA = oint x^{a+1} y^b/(a+1) dy on each edge, because the published (14f) does not reproduce the elementary case: for the CCW unit square it returns -1/3 where int y^2 dA = +1/3, i.e. its edge factor is (x_v - x_{v+1}) and not the printed (x_{v+1} - x_v). (14a)-(14e) do check out — (14d) agrees term for term with the derivation after the identity (x0+x1)(x0^2+x1^2) = x0^3 + x0^2 x1 + x0 x1^2 + x1^3. Every form below was verified against the unit square and, in the ctest, against a randomized triangulated oracle.

The whole PV assembly is quadratic in s (both A_ij = sum w s_i s_j and b_i = sum w s_i (b.s)), so a globally flipped polygon orientation cancels exactly; the sign convention below is nevertheless CCW-positive so s[0] is a readable area.

Definition at line 593 of file curvature.hpp.

Referenced by pvTermPolygon().

◆ ptFitInit()

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

Definition at line 643 of file curvature.hpp.

◆ ptFitAdd() [1/2]

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.

Parameters

X	the point in target-centred cell units.
org	the fit origin (the target cell's PLIC centroid), same units.
dW	Wendland support width in cell units.

Definition at line 656 of file curvature.hpp.

References peclet::core::vof::VofMetric::toPhys(), and wendlandWeight().

Referenced by ptFitAdd().

◆ ptFitAdd() [2/2]

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.

Definition at line 680 of file curvature.hpp.

References ptFitAdd().

◆ ptFitSolve()

bool peclet::core::vof::ptFitSolve	(	const PtFit &	f,
		double	a[6],
		bool &	red
	)

Solve the point fit, with the same rank cascade as pvFitSolve.

Definition at line 687 of file curvature.hpp.

References curvSolveSym().

◆ pvFitInit()

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

Definition at line 723 of file curvature.hpp.

◆ pvTermNormal()

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.

Returns false exactly where pvFitTerm returns false BEFORE it builds the polygon (np is then unspecified).

Definition at line 762 of file curvature.hpp.

References peclet::core::vof::VofMetric::h.

Referenced by pvFitTerm().

◆ pvTermPolygon()

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.

Sets t.ok and returns it. Requires nv >= 3.

pvFitTerm is pvTermNormal -> plicPolygon -> pvTermPolygon with the same arithmetic in the same order, so a caller that builds each cell's polygon ONCE (a per-cell cache, PvPolygon) and feeds it here reproduces pvFitTerm bit for bit (core tests/test_vof_pvcache.cpp T1, T2).

Definition at line 791 of file curvature.hpp.

References polygonMoments2d(), peclet::core::vof::VofMetric::toPhys(), and wendlandWeight().

Referenced by pvFitTerm().

◆ pvFitTerm()

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.

Same parameters as pvFitAdd. Sets t.ok and returns it: false exactly where pvFitAdd returns false (a rejected polygon), and then w, B, s are unspecified.

The composition pvTermNormal -> plicPolygon -> pvTermPolygon (design G §5.1a), bitwise the pre-split body (core tests/test_vof_pvcache.cpp T1 against a frozen verbatim copy).

Definition at line 847 of file curvature.hpp.

References plicPolygon(), pvTermNormal(), and pvTermPolygon().

Referenced by pvFitAdd().

◆ pvFitAccum()

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.

Calling it in the stencil's canonical order reproduces pvFitAdd's sums.

Definition at line 864 of file curvature.hpp.

Referenced by pvFitAdd().

◆ pvFitAdd() [1/2]

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.

Parameters as documented above PvTerm.

Returns: true if the polygon contributed.

Definition at line 878 of file curvature.hpp.

References pvFitAccum(), and pvFitTerm().

Referenced by pvFitAdd().

◆ pvFitAdd() [2/2]

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.

Definition at line 890 of file curvature.hpp.

References pvFitAdd().

◆ pvFitSolve()

bool peclet::core::vof::pvFitSolve	(	const PvFit &	f,
		double	a[6],
		bool &	red
	)

Solve the accumulated PV system.

A is a sum of npoly rank-1 outer products, so it cannot have rank 6 with fewer than 6 polygons and it is only generically full rank with more. On a rank-deficient system the fit drops to the 3-parameter sub-model ‘{1, x’^2, y'^2}— the exact minimiser of the SAME cost function restricted toa1 = a2 = a4 = 0`, which is the right restriction here because the frame is already aligned with the interface normal at the origin, so the linear terms are the small ones. Below three polygons there is nothing to fit and the caller must say so out loud.

Parameters

red	[out] true when the reduced model was used.

Returns: false when no fit was possible at all.

Definition at line 908 of file curvature.hpp.

References curvSolveSym().

◆ pvFitAccumLower()

void peclet::core::vof::pvFitAccumLower	(	PvFit &	f,
		const PvTerm &	t
	)

Lower-triangle accumulation: pvFitAccum restricted to j <= i (design G §5.1c).

The entries A[i][j], j > i, stay at pvFitInit's 0.0.

Why it is bitwise for the solve. curvSolveSym reads only the diagonal and the strict lower triangle of its input after the Jacobi scaling (it scales the upper half too, but nothing reads it), and pvFitSolve's reduced model copies A[3][0], A[5][0], A[5][3] — also lower — into the lower half of its 3x3 system. Every entry it does keep is the same sequential sum from 0.0 in the same order with the same expression, so pvFitSolve returns the same bits (core tests/test_vof_pvcache.cpp T3).

Definition at line 950 of file curvature.hpp.

◆ pvFitAccumEntry()

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.

The expressions are pvFitAccum's, so a per-entry loop over the terms in the same order equals pvFitAccum's interleaved loop bit for bit.

Definition at line 970 of file curvature.hpp.

◆ pvFitFromEntries()

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.

Definition at line 985 of file curvature.hpp.

◆ pvPolygonBuild()

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.

Definition at line 1004 of file curvature.hpp.

References plicPolygon().

◆ pvOutsideSupport()

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.

Why skipping it is bitwise. pvTermPolygon's radius is the norm of the frame projections of the mean of X_k = toPhys(off + v_k - 0.5 - org); this is the norm of the same vector, computed in a different order, so the two differ by a few ulp (~1e-15 relative). The margin (1 + 1e-9) on r^2 is 5e-10 on r, so every skipped polygon has r > dW, hence q >= 1 and w = 0: the accepted set, npoly and every accumulated byte are unchanged (core tests/test_vof_pvcache.cpp T2: false skips must be 0). A polygon with nv < 3 is never skipped (the polygon half rejects it on its own), nor is anything when dW <= 0, where wendlandWeight is 1 everywhere (an unbounded support).

Definition at line 1033 of file curvature.hpp.

References peclet::core::vof::VofMetric::toPhys().

◆ pvMomentsBuild()

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.

For the triangle (A, B, C) with area At and S = A + B + C: int dA = At, int Y dA = At S / 3, and int Y_p Y_q dA = At/12 (A_p A_q + B_p B_q + C_p C_q + S_p S_q) — the exact triangle moments (int l_p l_q dA = At (1 + d_pq) / 12 in barycentric coordinates). The polygon is convex, so every At >= 0.

Definition at line 1063 of file curvature.hpp.

References plicPolygon(), peclet::core::vof::VofMetric::toPhys(), and vofPhysNormalInv().

◆ pvFrameMoments()

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).

false when there is no normal. Against curvFallbackFrame: the frame vectors are bitwise (same vofPhysNormalInv, same curvFrame), and org is toPhys of its index-space origin up to rounding, because the area centroid maps linearly under H (T4: |d org| <= 1e-14).

Definition at line 1121 of file curvature.hpp.

References curvFrame(), and peclet::core::vof::VofMetric::toPhys().

◆ pvTermMoments()

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.

off is the stencil cell's integer index offset, org the target's PHYSICAL origin (pvFrameMoments). Sets t.ok and returns it.

Equivalence. With X = Y + d relative to the target origin (d = toPhys(off) - org), the projected polygon's moments are ‘int_{P’} phi(x', y') dA' = (n_j . n_t) int_P phi(t1.X, t2.X) dA: projecting a planar region scales its area by the cosine, and the projected polygon is the polygon of the projected vertices. Andint_{P'} z' dA' = (n_j . n_t) n_t . int_P X dA, because X lies on the plane — today'sB = b0 s0 + b1 s1 + b2 s2is that same integral. The Wendland point (the vertex average of X_k inpvFitTerm) isP.c + d. So the term differs from pvFitTerm's only by the rounding sequence (and by the orientation sign of the projected polygon, which cancels in the quadratic normal equations: heres0 = cj a > 0`). A RECORDED numerics change at round-off (max rel. d kappa 7.7e-15, design G §2.3).

Definition at line 1154 of file curvature.hpp.

References peclet::core::vof::VofMetric::toPhys(), and wendlandWeight().

◆ vofEpsEff()

double peclet::core::vof::vofEpsEff	(	double	eps,
		bool	solid
	)

The effective fluid volume of a cell (rule 1).

Returns 0 for a solid cell, so a caller can use eps_eff > 0 as the "is advected" predicate and never divide by zero.

Definition at line 101 of file cutcell.hpp.

References kVofEpsFloor.

◆ vofIsSolidCell()

bool peclet::core::vof::vofIsSolidCell	(	double	eps,
		const double	oLo[3],
		const double	oHi[3]
	)

Is this cell SOLID? eps == 0 and every one of the six faces closed.

oLo/oHi are the openness of the -d/+d faces.

Definition at line 109 of file cutcell.hpp.

◆ vofCutFlux()

double peclet::core::vof::vofCutFlux	(	double	openness,
		double	rawFlux
	)

The openness-weighted geometric flux (1): one number per face, used with opposite signs by the two neighbours, so the sum telescopes bit-exactly whatever o is.

Definition at line 120 of file cutcell.hpp.

◆ vofCutFluxClamp()

double peclet::core::vof::vofCutFluxClamp	(	double	fluxMag,
		double	sweep,
		double	epsDon,
		double	cDon
	)

Weymouth's admissible flux interval, generalized to a cut donor cell — the cut-cell twin of momentum_advect.hpp point 3, and for the same reason.

The whole-cell PLIC slab times the open area is an APPROXIMATION of the solid-clipped flux polygon, so it is not automatically bounded by what the donor actually holds: through a face of open area o and slab thickness |a| the scheme sweeps a FLUID volume o|a|, of which at most eps_don C_don can be liquid and at most eps_don (1 - C_don) can be gas. Hence

max(0, o|a| - eps_don (1 - C_don))  <=  |F|  <=  min(o|a|, eps_don C_don).

It is applied to the ONE face value both neighbours share, so conservation still telescopes bit-exactly; it is inactive wherever the geometry is self-consistent; and it is applied ONLY when the donor is MIXED — a pure-phase donor takes wyFaceFlux's algebraic branch, whose flux is o C_don |a| and therefore already exactly bounded, and clamping it would break the exact full-cell cancellation (o|a| vs min(o|a|, eps)) that makes an interior full cell stationary to the last bit.

MEASURED without it (24^3 packing, 200 kinematic steps at CFL 0.2): the [0,1] clip of rule 3 fires at up to 3.2e-5 liquid volume per step and the conserved functional drifts 1.3e-8 in 30 steps. With it: the clip stops firing and the drift is the projection floor.

Definition at line 144 of file cutcell.hpp.

◆ vofCutCourant()

double peclet::core::vof::vofCutCourant	(	double	openness,
		double	a,
		double	eps
	)

The boundedness-relevant Courant number of a face against a cell (rule 2).

Two constraints, and the max of them is the answer. WO-Q specifies o_f |a_f| / max(eps,0.1) — the flux volume against the cell's FLUID volume — and that is the second term. On its own it is WRONG, and loudly: it is smaller than |a_f| wherever o_f < eps, so a nearly-closed face (o = 0.02) inside an open cell would license |a| = 10. But |a| is the thickness of the slab the geometric flux clips out of the DONOR CELL, and a slab thicker than the cell is not a flux at all — plicSlabVolume is only a flux for |a| <= 1, and Weymouth's proof needs |a| <= 1/4 in 3D whatever the geometry is. MEASURED with the second term alone: a 24^3 sphere packing at "CFL 0.2" ran at dt = 1.85 and lost 70 % of the liquid volume in 200 steps while the flux sum still telescoped — the classic signature of an over-CFL Weymouth-Yue run (conservation is algebraic and survives; boundedness does not).

So: max( |a_f| , o_f |a_f| / max(eps_i, 0.1) ). It reduces EXACTLY to the uncut |a_f| in clear fluid (o = eps = 1), it throttles by 1/eps in a genuinely cut cell, and the 0.1 floor stops it chasing slivers to dt = 0 (a sliver's boundedness is handled by the clip, which is a measured number).

Definition at line 175 of file cutcell.hpp.

References kVofCourantEpsFloor.

◆ vofClipCut()

double peclet::core::vof::vofClipCut	(	double	c,
		double	epsEff,
		double &	moved
	)

The clip of rule 3, returned as the CLIPPED value; moved receives the liquid VOLUME the clip created (positive) or destroyed (negative) in this cell, i.e.

eps_eff * (C_clipped - C_raw).

Definition at line 184 of file cutcell.hpp.

◆ vofFillReadable()

bool peclet::core::vof::vofFillReadable	(	unsigned char	state,
		int	k
	)

Is state readable by pass k (1-based)? Fluid always; a solid cell only if an EARLIER pass filled it.

Definition at line 203 of file cutcell.hpp.

References kVofFillNone, and kVofFillPass0.

◆ vofGhostDepth1()

int peclet::core::vof::vofGhostDepth1	(	int	i,
		int	g,
		int	n
	)

Ghost depth of an index on one axis of an extended block: 0 for an inner cell, growing outward.

Definition at line 208 of file cutcell.hpp.

Referenced by vofGhostDepth().

◆ vofGhostDepth()

int peclet::core::vof::vofGhostDepth	(	int	x,
		int	y,
		int	z,
		int	g,
		int	nx,
		int	ny,
		int	nz
	)

Ghost depth of a cell = the max over the three axes.

Definition at line 217 of file cutcell.hpp.

References vofGhostDepth1().

◆ plicSq()

double peclet::core::vof::plicSq ( double x )

Definition at line 62 of file plic.hpp.

Referenced by plicAlphaReduced(), and plicVolumeReduced().

◆ plicCube()

double peclet::core::vof::plicCube ( double x )

Definition at line 65 of file plic.hpp.

Referenced by plicAlphaReduced(), and plicVolumeReduced().

◆ plicNormalizeL1()

double peclet::core::vof::plicNormalizeL1 ( double m[3] )

L1-normalize a normal in place. Returns the L1 norm before normalization (0 => untouched).

Definition at line 70 of file plic.hpp.

Referenced by mycNormal(), vofWettingPlane(), youngsNormal(), and youngsNormalFluidOnly().

◆ vofPhysNormal()

double peclet::core::vof::vofPhysNormal	(	const double	m[3],
		const VofMetric &	g,
		double	n[3]
	)

The PHYSICAL unit normal of an index-space normal m (any scale), i.e.

H^-1 m L2-normalized, and the shape factor s(m) = |H^-1 m| / |m| as the return value.

s is the ratio of two L2 norms of the SAME vector when H = I, so it is exactly 1.0 there and n is m/|m| — today's expression, bit for bit. Returns 0 (leaving n untouched) for a degenerate normal, exactly as pcUnitNormal does.

Why there are TWO of these, and it is not redundancy. x / q and x * (1/q) are DIFFERENT doubles in IEEE-754 (one rounding versus two), and the call sites this replaces are split between the two idioms: pcUnitNormal and the wetting frames divide, while the curvature and interface-area frames precompute invn = 1/|m| and multiply. Rule B (bit-identity on the isotropic path, flow/doc/anisotropic_vof.md §1) therefore requires each site to keep the form it had — so this one divides and vofPhysNormalInv multiplies by the reciprocal. Picking one for both would move the last bit of every curvature in the suite.

Definition at line 143 of file plic.hpp.

References peclet::core::vof::VofMetric::h.

Referenced by vofWettingPlane().

◆ vofPhysNormalInv()

double peclet::core::vof::vofPhysNormalInv	(	const double	m[3],
		const VofMetric &	g,
		double	n[3]
	)

vofPhysNormal in the RECIPROCAL-MULTIPLY form — n = (H^-1 m) * (1/|H^-1 m|), which is what the curvature cascade's and the interface area's frame construction has always computed.

See the note above for why the two forms both exist.

Definition at line 160 of file plic.hpp.

References peclet::core::vof::VofMetric::h.

Referenced by pvMomentsBuild().

◆ vofIndexNormal()

void peclet::core::vof::vofIndexNormal	(	const double	n[3],
		const VofMetric &	g,
		double	m[3]
	)

The INDEX-space normal of a physical normal n (any scale): m = H n, returned unnormalized (plicAlpha / plicVolume renormalize internally, so no L1 pass is needed).

Definition at line 177 of file plic.hpp.

References peclet::core::vof::VofMetric::h.

Referenced by vofWettingPlane().

◆ plicVolumeReduced()

double peclet::core::vof::plicVolumeReduced	(	double	w,
		double	n1,
		double	n2,
		double	n3
	)

Reduced-symmetry forward relation.

Inputs are already mirrored, sorted and L1-normalized: 0 <= n1 <= n2 <= n3, n1 + n2 + n3 = 1, and 0 <= w <= 1/2 is the plane offset measured from the origin corner. Returns the truncated volume in [0, 1/2].

Lehmann & Gekle (2022) eq. (11) / Listing 1. Case order (5) -> (1) -> (2) -> (3|4) is the paper's: it keeps the divisions by n1 (which may be zero for a 1D/2D normal) behind case conditions that are unreachable when n1 = 0.

Definition at line 194 of file plic.hpp.

References plicCube(), and plicSq().

Referenced by plicVolume().

◆ plicVolume()

double peclet::core::vof::plicVolume	(	double	mx,
		double	my,
		double	mz,
		double	alpha
	)

Volume fraction of the unit cell [0,1]^3 lying on the fluid side m .

x < alpha.

Parameters

mx,my,mz	interface normal (any scale; canonical form has \|mx\|+\|my\|+\|mz\| = 1)
alpha	plane offset in the same scale as m

Returns: V in [0,1]

Definition at line 221 of file plic.hpp.

References plicVolumeReduced().

Referenced by planeCellFraction(), planeCellFractionAniso(), plicBoxVolume(), plicSlabVolume(), and vofWettingFraction().

◆ plicAlphaReduced()

double peclet::core::vof::plicAlphaReduced	(	double	v,
		double	n1,
		double	n2,
		double	n3
	)

Reduced-symmetry analytic inverse of plicVolumeReduced.

0 <= n1 <= n2 <= n3, n1 + n2 + n3 = 1, 0 <= V <= 1/2; returns the offset w in [0, 1/2].

Lehmann & Gekle (2022) Listing 4 — the Scardovelli-Zaleski analytic solution in L1 form, branch-reduced. Case order (5) -> (2) -> (1) -> (3|4): every division by n2 or n1 sits behind a condition that already excluded the corresponding degenerate normal.

Definition at line 263 of file plic.hpp.

References plicCube(), and plicSq().

Referenced by plicAlpha().

◆ plicAlpha()

double peclet::core::vof::plicAlpha	(	double	mx,
		double	my,
		double	mz,
		double	v
	)

Plane offset alpha such that plicVolume(mx,my,mz,alpha) == V, exactly (analytic inverse).

Parameters

mx,my,mz	interface normal (any scale; alpha comes back in the same scale)
v	target volume fraction, clamped to [0,1]

Definition at line 303 of file plic.hpp.

References plicAlphaReduced().

Referenced by vofWettingPlane().

◆ plicSlabVolume()

double peclet::core::vof::plicSlabVolume	(	double	mx,
		double	my,
		double	mz,
		double	alpha,
		int	dir,
		double	a,
		double	b
	)

Fluid volume of the PLIC polyhedron inside the axis-aligned slab a <= x_dir <= b, expressed as a fraction of the WHOLE cell (so plicSlabVolume(..., 0, 1) == plicVolume(...)).

No clipping code: the slab is the unit cell of a rescaled coordinate x_dir = a + (b-a) xi, under which the plane becomes ‘m’ . y = alpha'withm'_dir = m_dir (b-a), the other components unchanged, andalpha' = alpha - m_dir a. The result is(b-a) * plicVolume(m')`.

Parameters

dir	0 = x, 1 = y, 2 = z

Definition at line 339 of file plic.hpp.

References plicVolume().

Referenced by faceFluxVolume().

◆ plicBoxVolume()

double peclet::core::vof::plicBoxVolume	(	double	mx,
		double	my,
		double	mz,
		double	alpha,
		double	a0,
		double	b0,
		double	a1,
		double	b1,
		double	a2,
		double	b2
	)

Fluid volume of the PLIC polyhedron inside the axis-aligned BOX [a0,b0] x [a1,b1] x [a2,b2], as a fraction of the WHOLE cell.

The three-axis generalization of plicSlabVolume and, like it, a pure coordinate rescale rather than a polyhedron clipper: under x_k = a_k + (b_k - a_k) xi_k the plane becomes ‘m’.xi = alpha'withm'_k = m_k (b_k - a_k) andalpha' = alpha - sum_k m_k a_k, so the answer isprod_k (b_k - a_k) * plicVolume(m')`.

Added at rung V2b (WO-K): the momentum control volume of component e is shifted half a cell in e, so its flux slab for a sweep in a TRANSVERSE direction d is a box that is half a cell wide in e and |a| wide in d — two axes at once, which plicSlabVolume cannot express. plicBoxVolume(..., 0,1, 0,1, 0,1) == plicVolume(...) bitwise, and the single-axis restriction reproduces plicSlabVolume bitwise (both gated in tests/kokkos/test_vof_plic.cpp).

Definition at line 361 of file plic.hpp.

References plicVolume().

◆ faceFluxVolume()

double peclet::core::vof::faceFluxVolume	(	double	mx,
		double	my,
		double	mz,
		double	alpha,
		int	dir,
		double	f
	)

Donor-cell geometric flux volume: the fluid volume swept out of the slab [0, f] along dir, as a fraction of the whole cell.

f = 1 reproduces plicVolume bitwise.

Definition at line 374 of file plic.hpp.

References plicSlabVolume().

◆ plicSt()

int peclet::core::vof::plicSt	(	int	i,
		int	j,
		int	k
	)

Stencil index: x-fastest, offsets (i-1, j-1, k-1) with i,j,k in {0,1,2} (suite/docs/CONVENTIONS.md).

c[plicSt(1,1,1)] is the centre cell.

Definition at line 385 of file plic.hpp.

◆ youngsNormal()

void peclet::core::vof::youngsNormal	(	const double	c[27],
		double	m[3]
	)

Youngs' normal: the 27-point weighted gradient (corner 1, edge 2, face 4 on each of the two opposing planes), L1-normalized, pointing into the gas (m ~ -grad C).

Degenerate (zero gradient) stencils return (1,0,0).

Definition at line 394 of file plic.hpp.

References PECLET_VOF_C, and plicNormalizeL1().

Referenced by mycNormal().

◆ mycNormal()

void peclet::core::vof::mycNormal	(	const double	c[27],
		double	m[3]
	)

MYC — Mixed Youngs-Centred normal (Aulisa et al.

2007). Three "centred column" candidates (one per axis: the plane written as sgn(m_d) x_d = ... + alpha, with the transverse components taken from central differences of the 3-cell column sums) plus the Youngs candidate. The centred candidate with the largest dominant component |m[d][d]| is selected first; Youngs replaces it when that dominant component exceeds Youngs' own largest component — the signature of a saturated (clipped) column, where the height difference under-reports the true slope and the centred estimate is artificially flat.

Transcribed from basilisk/src/myc.h. The returned normal is L1-normalized and points into the gas. As noted there, the centre cell never enters either estimator, so an isolated droplet gives a degenerate stencil; that case returns (1,0,0).

Definition at line 451 of file plic.hpp.

References PECLET_VOF_C, plicNormalizeL1(), and youngsNormal().

◆ planeCellFraction()

double peclet::core::vof::planeCellFraction	(	double	mx,
		double	my,
		double	mz,
		double	alphaGlobal,
		double	x0,
		double	y0,
		double	z0,
		double	h
	)

Exact fluid fraction of the axis-aligned cell [x0,x0+h] x [y0,y0+h] x [z0,z0+h] under the GLOBAL plane m .

X = alphaGlobal (fluid on the m . X < alphaGlobal side). This is just plicVolume on the cell-local plane, so initPlane data is exact to round-off.

Definition at line 542 of file plic.hpp.

References plicVolume().

Referenced by sphereCellFraction().

◆ planeCellFractionAniso()

double peclet::core::vof::planeCellFractionAniso	(	double	mx,
		double	my,
		double	mz,
		double	alphaGlobal,
		double	x0,
		double	y0,
		double	z0,
		const VofMetric &	g
	)

The same for an anisotropic cell [x0, x0+g.h[0]] x [y0, y0+g.h[1]] x [z0, z0+g.h[2]]: the index normal of the global plane is H n (vofIndexNormal), so the exact fraction is one plicVolume call.

Reproduces planeCellFraction(..., h) bitwise for g.h = {h,h,h} (the same three products in the same order) and is the oracle of the K1 gate.

Definition at line 552 of file plic.hpp.

References peclet::core::vof::VofMetric::h, and plicVolume().

◆ sphereCellFraction()

double peclet::core::vof::sphereCellFraction	(	double	cx,
		double	cy,
		double	cz,
		double	r,
		double	x0,
		double	y0,
		double	z0,
		double	h,
		int	levels
	)

Fluid fraction of a cell for a sphere (fluid = INSIDE the sphere), by recursive octree subdivision with an exact tangent-plane closure at the leaves.

Sub-boxes entirely inside or entirely outside short-circuit on the centre-distance / half-diagonal test, so only the O(4^L) boxes the surface actually crosses are refined. Iterative (explicit stack) — device-safe, no recursion.

The leaf closure is second order in the leaf size delta, giving a cell-fraction error ~ delta^2 / (8 h R): levels = 4 is already ~3e-5 at R/h = 16. Test-only helper.

Parameters

levels subdivision depth, clamped to [0, 8]

Definition at line 569 of file plic.hpp.

References planeCellFraction().

◆ youngsNormalFluidOnly()

bool peclet::core::vof::youngsNormalFluidOnly	(	const double	c[27],
		const unsigned char	fluid[27],
		double	m[3]
	)

Youngs' 27-point weighted gradient restricted to non-solid cells.

Youngs' weights on each of the two opposing planes of axis d are 4 (plane centre), 2 (edge) and 1 (corner), i.e. w = (2 - |t1|)(2 - |t2|) over the two transverse offsets, and sum to 16. Here each half-plane sum is taken over the INCLUDED (fluid) cells only and divided by the included weight, so the estimator is a difference of two weighted MEANS and the omission of a solid cell does not read as a colour deficit. (With every cell fluid this is Youngs' own normal divided by 16 — the same direction, hence the same PLIC plane.)

If one of the two planes has no fluid cell at all — the normal case for the wall-normal axis right next to a wall — the difference is taken one-sided against the CENTRE plane and doubled, which is the same first-order slope over half the distance. If both outer planes are empty the component is 0.

Parameters

c	27-cell colour stencil, `plicSt` order
fluid	1 where the cell is fluid (not `kVofSolid`), 0 where it is solid
m	out: L1-normalized normal pointing into the gas

Returns: false if the gradient is degenerate (no usable direction) — the caller must then fall back to the neutral rule rather than invent one.

Definition at line 166 of file wetting.hpp.

References plicNormalizeL1().

◆ vofPlicCentroid()

void peclet::core::vof::vofPlicCentroid	(	const double	m[3],
		double	alpha,
		double	p[3]
	)

The centroid of the PLIC polygon of cell f in the cell's own unit-cube frame, for the plane m .

x = alpha. Falls back to the plane point nearest the cell centre when the polygon is degenerate (fewer than three vertices), which is the WO-S fallback written out.

Definition at line 195 of file wetting.hpp.

References plicPolygon(), and polygonAreaCentroid().

Referenced by vofWettingPlane().

◆ vofWettingPlane() [1/2]

int peclet::core::vof::vofWettingPlane	(	const double	mfIn[3],
		double	cf,
		const double	nwIn[3],
		double	cosT,
		double	sinT,
		double	sdfF,
		int	pivot,
		double	tEps,
		double	mth[3],
		double &	alphaTh,
		double &	cosApp,
		const VofMetric &	g
	)

Build the theta-plane in the FLUID cell's own unit-cube frame.

Parameters

mfIn	the fluid-side normal (any scale; `youngsNormalFluidOnly`'s output). Only its component IN the wall is used by the default anchor — see the file header.
cf	the fluid cell's colour (reconstructed on the WHOLE unit cell, as WO-Q's flux is)
nwIn	the wall normal, solid -> fluid (any scale)
cosT	cos(theta),
sinT	sin(theta) (theta prescribed, through the liquid)
sdfF	the signed distance at the fluid cell's CENTRE, in cells (only read by the `kVofPivotWallNormal` / `kVofPivotContactLine` ablations)
pivot	`VofWettingPivot`
tEps	below `\|t\| = \|sin(theta_apparent)\| = tEps` no rotation is defined
mth	out: the theta-plane normal, L2-normalized, pointing into the gas
alphaTh	out: its offset in the fluid cell's frame (liquid side `mth . x < alphaTh`)
cosApp	out: the measured apparent cosine `m_f . n_w` (the diagnostic G1 reports)

Returns: the branch (kVofWetTheta or kVofWetParallel); never fails.

Definition at line 226 of file wetting.hpp.

References peclet::core::vof::VofMetric::h, kVofPivotContactLine, kVofPivotVolume, kVofPivotWallNormal, kVofWetParallel, kVofWetTheta, plicAlpha(), plicNormalizeL1(), peclet::core::vof::VofMetric::toPhys(), vofIndexNormal(), vofPhysNormal(), and vofPlicCentroid().

Referenced by vofWettingPlane().

◆ vofWettingPlane() [2/2]

int peclet::core::vof::vofWettingPlane	(	const double	mfIn[3],
		double	cf,
		const double	nwIn[3],
		double	cosT,
		double	sinT,
		double	sdfF,
		int	pivot,
		double	tEps,
		double	mth[3],
		double &	alphaTh,
		double &	cosApp
	)

Unit-metric overload (isotropic cells) — the pre-Phase-3 signature, unchanged arithmetic.

Definition at line 314 of file wetting.hpp.

References vofWettingPlane().

◆ vofWettingFraction()

double peclet::core::vof::vofWettingFraction	(	const double	mth[3],
		double	alphaTh,
		const int	ds[3]
	)

The liquid fraction of the solid cell at integer offset ds from the fluid cell, under the theta-plane returned by vofWettingPlane.

Clamped to [0,1].

Definition at line 324 of file wetting.hpp.

References plicVolume().

Variable Documentation

◆ kHfColumn

constexpr int peclet::core::vof::kHfColumn = 7

inlineconstexpr

Column length of the height function, NH in Han et al.

eq. (3). 7 is the value used throughout the literature (Popinet 2009; Cummins 2005; Owkes 2015; Jibben 2019) and is exactly what the colour field's g = 3 halo reaches.

Definition at line 148 of file curvature.hpp.

◆ kPvHalf

constexpr int peclet::core::vof::kPvHalf = 2

inlineconstexpr

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

Definition at line 150 of file curvature.hpp.

◆ kPvWeightWidth

constexpr double peclet::core::vof::kPvWeightWidth = 2.5

inlineconstexpr

Wendland support width d in cell units.

Han et al. §5 use d = 2.5 with S = 5 and show that d = 3.5 recovers first-order convergence of the spurious currents on a translating droplet while d = 4.5 destroys convergence entirely (over-smoothing). Exposed as an argument so the V4 rung can sweep it; 2.5 is the default and the static-test value.

Definition at line 155 of file curvature.hpp.

◆ kPtWeightWidth

constexpr double peclet::core::vof::kPtWeightWidth = 3.5

inlineconstexpr

Wendland support width for the tier-2b height-position fit.

Wider than the PV one because the points it fits reach sqrt(1 + 1 + 2.5^2) ~ 2.9 cells by construction and a 2.5-cell support would silently discard the very columns that made the tier fire.

Definition at line 159 of file curvature.hpp.

◆ kHfPureEps

constexpr double peclet::core::vof::kHfPureEps = 1e-10

inlineconstexpr

The FLOOR of the "this cell is pure" thresholds.

A cell at 1 - 1e-12 is full for the purposes of a column-consistency test; the residue it contributes to the column sum is far below the height's own discretisation error.

Definition at line 164 of file curvature.hpp.

Referenced by hfIsEmpty(), and hfIsFull().

◆ kPvEntries

constexpr int peclet::core::vof::kPvEntries = 27

inlineconstexpr

Number of independent entries of the PV normal equations: the 21 of the lower triangle (with the diagonal) and the 6 of b.

Definition at line 963 of file curvature.hpp.

◆ kVofEpsFloor

constexpr double peclet::core::vof::kVofEpsFloor = 1.0 / 64.0

inlineconstexpr

The smallest fluid volume buildCellFraction's 4³ subsampling can resolve (one subcell of 64).

Used as the floor of the effective fluid volume the cut-cell update divides by — see rule 1.

Definition at line 85 of file cutcell.hpp.

Referenced by vofEpsEff().

◆ kVofCourantEpsFloor

constexpr double peclet::core::vof::kVofCourantEpsFloor = 0.1

inlineconstexpr

Below this the interface-local Courant number stops being amplified by 1/eps (rule 2): a cut cell at eps = 0.1 already reports 10x its face Courant number, and chasing thinner slivers with dt buys boundedness the clip provides for free.

Definition at line 90 of file cutcell.hpp.

Referenced by vofCutCourant().

Classes

Enumerations

Functions

Variables

Enumeration Type Documentation

◆ CurvatureBranch

◆ VofCellKind

◆ VofFillState

◆ VofWettingPivot

◆ VofWettingBranch

Function Documentation

◆ hfIsFull()

◆ hfIsEmpty()

◆ hfColumnHeight()

◆ hfPatchKappa() [1/2]

◆ hfPatchKappa() [2/2]

◆ curvFrame()

◆ wendlandWeight()

◆ paraboloidKappa()

◆ curvSolveSym()

◆ plicPolygon()

◆ polygonAreaCentroid()

◆ polygonMoments2d()

◆ ptFitInit()

◆ ptFitAdd() [1/2]

◆ ptFitAdd() [2/2]

◆ ptFitSolve()

◆ pvFitInit()

◆ pvTermNormal()

◆ pvTermPolygon()

◆ pvFitTerm()

◆ pvFitAccum()

◆ pvFitAdd() [1/2]

◆ pvFitAdd() [2/2]

◆ pvFitSolve()

◆ pvFitAccumLower()

◆ pvFitAccumEntry()

◆ pvFitFromEntries()

◆ pvPolygonBuild()

◆ pvOutsideSupport()

◆ pvMomentsBuild()

◆ pvFrameMoments()

◆ pvTermMoments()

◆ vofEpsEff()

◆ vofIsSolidCell()

◆ vofCutFlux()

◆ vofCutFluxClamp()

◆ vofCutCourant()

◆ vofClipCut()

◆ vofFillReadable()

◆ vofGhostDepth1()

◆ vofGhostDepth()

◆ plicSq()

◆ plicCube()

◆ plicNormalizeL1()

◆ vofPhysNormal()

◆ vofPhysNormalInv()

◆ vofIndexNormal()

◆ plicVolumeReduced()

◆ plicVolume()

◆ plicAlphaReduced()

◆ plicAlpha()

◆ plicSlabVolume()

◆ plicBoxVolume()

◆ faceFluxVolume()

◆ plicSt()

◆ youngsNormal()

◆ mycNormal()

◆ planeCellFraction()

◆ planeCellFractionAniso()

◆ sphereCellFraction()

◆ youngsNormalFluidOnly()

◆ vofPlicCentroid()

◆ vofWettingPlane() [1/2]

◆ vofWettingPlane() [2/2]

◆ vofWettingFraction()

Variable Documentation

◆ kHfColumn

◆ kPvHalf

◆ kPvWeightWidth