marching_cubes.hpp File Reference

flow — WO-P3d: a JOINED interfacial area, from marching tetrahedra on the cell-centre lattice. More...

#include <Kokkos_Core.hpp>
#include <Kokkos_MathematicalFunctions.hpp>

Go to the source code of this file.

Classes
struct	peclet::flow::vof::McVertex
	One dual-cube corner = one cell centre. More...

Namespaces
namespace	peclet
namespace	peclet::flow
namespace	peclet::flow::vof

Enumerations
enum	peclet::flow::vof::McSource : int { peclet::flow::vof::kMcSrcColour = 0 , peclet::flow::vof::kMcSrcPlic = 1 }
	What the edge crossing interpolates. More...
enum	peclet::flow::vof::McDeposit : int { peclet::flow::vof::kMcDepositCentroid = 0 , peclet::flow::vof::kMcDepositSplit = 1 }
	Which cell a triangle's area is booked to. More...

Functions
void	peclet::flow::vof::mcCornerPos (int k, double p[3])
	Local coordinates of dual-cube corner k (bit 0 = x, bit 1 = y, bit 2 = z).
double	peclet::flow::vof::mcEdgeT (const McVertex &a, const McVertex &b, const double dv[3], int src)
	The parameter t in [0,1] at which the level set cuts the edge from corner a to corner b.
double	peclet::flow::vof::mcTwiceArea (const double p0[3], const double p1[3], const double p2[3], const VofMetric &g)
	Twice the area of the triangle (p0, p1, p2) — i.e. |(p1-p0) x (p2-p0)|.
double	peclet::flow::vof::mcPolygonArea (const double p[8][3], int n, const VofMetric &g)
	Area of a PLANAR polygon of n vertices, 1/2 |sum p_i x p_{i+1}|.
int	peclet::flow::vof::mcClipHalf (double p[8][3], int n, int axis, double c, bool keepAbove)
	One Sutherland-Hodgman clip of a convex polygon against the half-space keepAbove ? x_axis >= c : x_axis <= c.
double	peclet::flow::vof::mcClipToCorner (const double p0[3], const double p1[3], const double p2[3], int k, const VofMetric &g)
	The area of the part of one triangle inside corner k's octant of the dual cube.
int	peclet::flow::vof::mcNearestMasked (const double p[3], int mask)
	The corner in mask nearest the point p; the lowest index breaks a tie.
int	peclet::flow::vof::mcRetarget (int k, int mask)
	The retarget rule.
double	peclet::flow::vof::mcTriangleToCorner (const double p0[3], const double p1[3], const double p2[3], int lc, int deposit, int mask, const VofMetric &g)
	The part of one triangle that belongs to local corner lc of the dual cube.
double	peclet::flow::vof::mcCubeCornerArea (const McVertex v[8], int lc, int src, int deposit, const VofMetric &g)
	The interfacial area inside one dual cube that belongs to its local corner lc.
double	peclet::flow::vof::mcCubeCornerArea (const McVertex v[8], int lc, int src, int deposit)
	Unit-metric overload (isotropic cells) — the pre-Phase-3 signature, unchanged arithmetic.

Detailed Description

flow — WO-P3d: a JOINED interfacial area, from marching tetrahedra on the cell-centre lattice.

Container-free by contract, exactly like plic.hpp, interface_area.hpp and phase_change.hpp: scalars and small local arrays only. The block-walking driver is vof/marching_cubes_field.hpp.

Why this file exists (the measurement, WO-P3c)

Every PER-CELL interfacial area — the PLIC polygon (kAreaPlic), the height function's metric (kAreaMetric), the plane rebuilt on the cascade normal (kAreaNormal), the height function's own footprint (kAreaFootprint) — is first order in h/R on a curved interface, and WO-P3c proved it with two analytic controls that contain none of the code under test: an exact-fraction circle with the EXACT radial normal and the analytic chord sums to −2.86 / −1.47 / −1.42 / −0.84 % of 2 pi R at R = 8/12/20/28, and the mode-3 footprint construction with EXACT heights sums to −3.50 / −2.42 / −1.66 / −1.11 / −0.81 % at R = 8/12/20/28/40, observed order 1.03. The mechanism is not the normal, the fractions, the metric or the footprint: a cell chooses its piece of the surface INDEPENDENTLY of its neighbours', so the pieces overlap and gap at every cell face instead of JOINING, and the mismatch is O(h^2 kappa) per face over O(L/h) faces, i.e. a relative O(h/R).

The cure is a surface whose pieces join by construction. That is what an isosurface triangulation is: one vertex per crossed lattice edge, shared by every cell that edge touches, so the triangles form a closed watertight sheet and the sum of their areas is the area of ONE surface rather than the sum of 6 % overlapping fragments. Marching cubes on the same colour fields was already measured 3-10x closer (WO-P3c: −0.59 / −0.49 % on the cylinder rows where every per-cell construction reads −1.4 … −1.6 %).

What is built here

Marching TETRAHEDRA on the Kuhn 6-tet decomposition of the dual cube (the 8 corners of a dual cube are 8 CELL CENTRES). Two reasons over the 256-case cube table:

• it is branch-light and table-free — 4 corner signs give one of three topologies (nothing, one triangle, one quad) and the whole case analysis is popcount — which is what makes it a clean device kernel, and
• the Kuhn decomposition is translation invariant, so two neighbouring dual cubes split their shared face on the SAME diagonal and the sheet is watertight with no ambiguous-face rule at all (the classic MC hole).

The two things that are chosen by MEASUREMENT, not taste

**(1) What is interpolated along an edge** (src):

• kMcSrcColour — the raw C = 1/2 level set, t = psi_a/(psi_a - psi_b) with psi = 1/2 - C. This is the work order's design, and it is REFUTED on a curved interface: +4.2 … +5.8 % of 4 pi R^2 over R = 4…28 with no convergence. The mechanism is the combination with the tetrahedra: C(d) is the SZ piecewise CUBIC of the centre distance, so linear interpolation of it is only accurate over a step short compared with a cell, and Kuhn's tets interpolate along sqrt(2) face diagonals and a sqrt(3) body diagonal. Every vertex on a long edge is misplaced along the normal, the sheet wrinkles, and a wrinkle only ever ADDS area — hence a one-signed 5 % that does not converge. (Marching CUBES, which interpolates only the 12 unit-length cube edges, reads −0.2 … +0.4 % on the same fields — the probe prints it as the external cross-check. So the defect is the combination, not either half alone.) It ships as the measured ablation.

• kMcSrcPlic — the zero of the PLIC-reconstructed signed distance. On ANY plane, tilted or not, both endpoint cells' planes ARE that plane, so every vertex lands exactly on it and the summed area is exact to round-off. This is the work order's named follow-on ("interpolate the PLIC-reconstructed signed distance instead of C") and it is what the 0.1 % tilted-plane gate needs.

  How the two endpoint planes are combined is itself a measurement, and the two candidates are three decades apart on the sphere probe. Averaging the two ROOTS reads +0.50 / +0.42 / +0.14 / +0.08 % of 4 pi R^2 at R = 8/12/20/28; blending the two signed distance FUNCTIONS, Phi(s) = (1-s) phi_a(s) + s phi_b(s) (a quadratic in s, solved exactly), reads +0.011 / +0.011 / +0.009 / +0.008 %. Both are exact on a plane; the blend wins because it weights each cell's plane by how NEAR the crossing is to it, so the Kuhn tets' long edges (sqrt(2) face diagonals, sqrt(3) the body diagonal) stop extrapolating a tangent plane most of a cell away from where it was fitted. The blend is what ships.

**(2) Which cell a piece of the sheet belongs to** (deposit):

• kMcDepositCentroid — the whole triangle goes to the cell containing its centroid.
• kMcDepositSplit — the triangle is CLIPPED to each cell's own cube and each cell gets the area of its piece. Both give the same SUM (the 8 cells of a dual cube tile it exactly), so they are indistinguishable on the a-priori area probe and differ only in the per-cell distribution the phase-change flux integral sees.

The deposit is a GATHER, never an atomic scatter

WO-P01's lesson: a bitwise-MPI phase change needs a fixed summation order. A triangle produced inside a dual cube can land in any of that cube's 8 corner cells, which as a SCATTER means atomics and a non-deterministic order. So the driver inverts it: a cell walks the 8 dual cubes it is a corner of, re-derives every triangle of each, and keeps only what belongs to itself. The 8x redundancy is the price of determinism, and it is what makes np 1/2/4 bitwise by construction rather than by tolerance.

Definition in file marching_cubes.hpp.