plic.hpp

Go to the documentation of this file.

 
#ifndef PECLET_CORE_VOF_PLIC_HPP
#define PECLET_CORE_VOF_PLIC_HPP
 
#include <Kokkos_Core.hpp>
#include <Kokkos_MathematicalFunctions.hpp>
 
namespace peclet::core::vof {
 
// ---------------------------------------------------------------------------------------------
// small local helpers
// ---------------------------------------------------------------------------------------------
 
KOKKOS_INLINE_FUNCTION double plicSq(double x) {
  return x * x;
}
KOKKOS_INLINE_FUNCTION double plicCube(double x) {
  return x * x * x;
}
 
KOKKOS_INLINE_FUNCTION double plicNormalizeL1(double m[3]) {
  const double s = Kokkos::fabs(m[0]) + Kokkos::fabs(m[1]) + Kokkos::fabs(m[2]);
  if (s > 0.0) {
    const double inv = 1.0 / s;
    m[0] *= inv;
    m[1] *= inv;
    m[2] *= inv;
  }
  return s;
}
 
// ---------------------------------------------------------------------------------------------
// the anisotropic cell metric (Phase 3 of suite/docs/PHYSICAL_UNITS_PLAN.md;
// flow/doc/anisotropic_vof.md §2)
// ---------------------------------------------------------------------------------------------
 
struct VofMetric {
  double h[3] = {1.0, 1.0, 1.0};
 
  KOKKOS_INLINE_FUNCTION bool isotropic() const { return h[0] == h[1] && h[1] == h[2]; }
  KOKKOS_INLINE_FUNCTION double det() const { return h[0] * h[1] * h[2]; }
  KOKKOS_INLINE_FUNCTION double maxH() const {
    const double m = h[0] > h[1] ? h[0] : h[1];
    return m > h[2] ? m : h[2];
  }
  KOKKOS_INLINE_FUNCTION double minH() const {
    const double m = h[0] < h[1] ? h[0] : h[1];
    return m < h[2] ? m : h[2];
  }
  KOKKOS_INLINE_FUNCTION void toPhys(const double X[3], double out[3]) const {
    out[0] = X[0] * h[0];
    out[1] = X[1] * h[1];
    out[2] = X[2] * h[2];
  }
};
 
KOKKOS_INLINE_FUNCTION double vofPhysNormal(const double m[3], const VofMetric& g, double n[3]) {
  const double q = Kokkos::sqrt(m[0] * m[0] + m[1] * m[1] + m[2] * m[2]);
  if (!(q > 0.0))
    return 0.0;
  const double p[3] = {m[0] / g.h[0], m[1] / g.h[1], m[2] / g.h[2]};
  const double pn = Kokkos::sqrt(p[0] * p[0] + p[1] * p[1] + p[2] * p[2]);
  if (!(pn > 0.0))
    return 0.0;
  n[0] = p[0] / pn;
  n[1] = p[1] / pn;
  n[2] = p[2] / pn;
  return pn / q;
}
 
KOKKOS_INLINE_FUNCTION double vofPhysNormalInv(const double m[3], const VofMetric& g, double n[3]) {
  const double q = Kokkos::sqrt(m[0] * m[0] + m[1] * m[1] + m[2] * m[2]);
  if (!(q > 0.0))
    return 0.0;
  const double p[3] = {m[0] / g.h[0], m[1] / g.h[1], m[2] / g.h[2]};
  const double n2 = p[0] * p[0] + p[1] * p[1] + p[2] * p[2];
  if (!(n2 > 0.0))
    return 0.0;
  const double invn = 1.0 / Kokkos::sqrt(n2);
  n[0] = p[0] * invn;
  n[1] = p[1] * invn;
  n[2] = p[2] * invn;
  return Kokkos::sqrt(n2) / q;
}
 
KOKKOS_INLINE_FUNCTION void vofIndexNormal(const double n[3], const VofMetric& g, double m[3]) {
  m[0] = n[0] * g.h[0];
  m[1] = n[1] * g.h[1];
  m[2] = n[2] * g.h[2];
}
 
// ---------------------------------------------------------------------------------------------
// (1) forward: plane -> volume
// ---------------------------------------------------------------------------------------------
 
KOKKOS_INLINE_FUNCTION double plicVolumeReduced(double w, double n1, double n2, double n3) {
  const double n12 = n1 + n2;
  // case (5) — the slab regime (the plane exits through both faces normal to axis 3), no division
  // by n1/n2. NOTE the condition is `n1+n2 <= w`, NOT Lehmann & Gekle's Listing 1 shortcut
  // `min(n1+n2,n3) <= w`: eq. (11) defines case (5) as "the remaining free sector mutually excluded
  // by the other four cases", which is the interval [n1+n2, n3] and is empty unless n3 >= n1+n2.
  // Hoisting the test to the front (as the listing does, for speed) makes `min(n1+n2,n3) <= w`
  // additionally fire at the single point w == n3 when n3 < n1+n2, where the correct case is (3):
  // e.g. m = (1/3,1/3,1/3), alpha = 1/3 gives 0 instead of the exact 1/6. Using the true case-(5)
  // interval keeps the front position (and its n1 = 0 protection) without the boundary defect.
  if (n12 <= w && w <= n3)
    return (w - 0.5 * n12) / n3;
  if (w < n1)  // case (1) — corner tetrahedron (implies n1 > 0)
    return plicCube(w) / (6.0 * n1 * n2 * n3);
  if (w <= n2)  // case (2) — one corner clipped; the n1 factor cancels, so n1 = 0 is safe
    return (3.0 * w * (w - n1) + plicSq(n1)) / (6.0 * n2 * n3);
  // cases (3) and (4) — two or three corners clipped; fdim() merges the two branches. Reachable
  // only with n1 > 0 (n1 = 0 is fully covered by cases (2) and (5) above).
  const double t3 = Kokkos::fmax(w - n3, 0.0);
  return (plicCube(w) - plicCube(w - n1) - plicCube(w - n2) - plicCube(t3)) / (6.0 * n1 * n2 * n3);
}
 
KOKKOS_INLINE_FUNCTION double plicVolume(double mx, double my, double mz, double alpha) {
  const double ax = Kokkos::fabs(mx), ay = Kokkos::fabs(my), az = Kokkos::fabs(mz);
  const double s = ax + ay + az;
  if (!(s > 0.0))
    return alpha > 0.0 ? 1.0 : 0.0;  // degenerate normal: the half-space is trivial
 
  // Mirror the negative components into the positive octant. Substituting x_i -> 1 - x_i turns
  // m_i x_i into m_i - m_i x_i', so the offset picks up |m_i| for every negative component.
  double beta = alpha;
  if (mx < 0.0)
    beta += ax;
  if (my < 0.0)
    beta += ay;
  if (mz < 0.0)
    beta += az;
 
  const double a1 = beta / s;  // offset in the L1-normalized, positive-octant frame
  if (a1 <= 0.0)
    return 0.0;
  if (a1 >= 1.0)
    return 1.0;
 
  // Sort ascending. n2 from the L1 identity (never from a third min/max) so n1+n2+n3 = 1 holds.
  const double n1 = Kokkos::fmin(Kokkos::fmin(ax, ay), az) / s;
  const double n3 = Kokkos::fmax(Kokkos::fmax(ax, ay), az) / s;
  const double n2 = Kokkos::fmax(1.0 - n1 - n3, 0.0);
 
  // V(alpha) + V(1-alpha) = 1 (point reflection x -> 1-x), so only w = min(a1, 1-a1) is computed.
  const double v = plicVolumeReduced(Kokkos::fmin(a1, 1.0 - a1), n1, n2, n3);
  return (a1 <= 0.5) ? v : 1.0 - v;
}
 
// ---------------------------------------------------------------------------------------------
// (2) inverse: volume -> plane
// ---------------------------------------------------------------------------------------------
 
KOKKOS_INLINE_FUNCTION double plicAlphaReduced(double v, double n1, double n2, double n3) {
  const double n12 = n1 + n2, n3v = n3 * v;
  if (n12 <= 2.0 * n3v)
    return n3v + 0.5 * n12;  // case (5); also catches n1 = n2 = 0 (1D normal)
 
  // past case (5), n2 > 0 (n2 = 0 would force n1 = 0 and hence n12 = 0 <= 2 n3 V)
  const double sqn1 = plicSq(n1), n26 = 6.0 * n2, v1 = sqn1 / n26;
  if (v1 <= n3v && n3v < v1 + 0.5 * (n2 - n1))  // case (2); also catches n1 = 0 (2D normal)
    return 0.5 * (n1 + Kokkos::sqrt(sqn1 + 8.0 * n2 * (n3v - v1)));
 
  const double v6 = n1 * n26 * n3v;  // = 6 n1 n2 n3 V
  if (n3v < v1)
    return Kokkos::cbrt(v6);  // case (1)
 
  // past case (2), n1 > 0. v3 is n3 * V at the (3)|(4) boundary w = n3 (or at the case (5)
  // boundary w = n1+n2 when the plane can never reach w = n3).
  const double v3 = (n3 < n12) ? (plicSq(n3) * (3.0 * n12 - n3) + sqn1 * (n1 - 3.0 * n3) +
                                  plicSq(n2) * (n2 - 3.0 * n3)) /
                                     (n1 * n26)
                               : 0.5 * n12;
  const double sqn12 = sqn1 + plicSq(n2), v6cb = v6 - plicCube(n1) - plicCube(n2);
  const bool case3 = (n3v < v3);
  // Both cubic cases reduce to  w^3 - 3 c w^2 + 3 b w + a = 0 :
  const double a = case3 ? v6cb : 0.5 * (v6cb - plicCube(n3));
  const double b = case3 ? sqn12 : 0.5 * (sqn12 + plicSq(n3));
  const double c = case3 ? n12 : 0.5;
  // c^2 - b > 0 is guaranteed here: 2 n1 n2 > 0 in case (3); 1/4 - (n1^2+n2^2+n3^2)/2 > 0 in case
  // (4), which is reachable only when n3 < n1+n2, i.e. n3 < 1/2.
  const double t = Kokkos::sqrt(c * c - b);
  // Casus irreducibilis: three real roots, |arg| <= 1 analytically. The clamp is a floating-point
  // guard on the case boundaries only (an out-of-range arg would otherwise return NaN).
  const double arg =
      Kokkos::fmin(Kokkos::fmax((plicCube(c) - 0.5 * a - 1.5 * b * c) / plicCube(t), -1.0), 1.0);
  return c - 2.0 * t * Kokkos::sin((1.0 / 3.0) * Kokkos::asin(arg));
}
 
KOKKOS_INLINE_FUNCTION double plicAlpha(double mx, double my, double mz, double v) {
  const double ax = Kokkos::fabs(mx), ay = Kokkos::fabs(my), az = Kokkos::fabs(mz);
  const double s = ax + ay + az;
  if (!(s > 0.0))
    return 0.0;
 
  const double n1 = Kokkos::fmin(Kokkos::fmin(ax, ay), az) / s;
  const double n3 = Kokkos::fmax(Kokkos::fmax(ax, ay), az) / s;
  const double n2 = Kokkos::fmax(1.0 - n1 - n3, 0.0);
 
  const double vc = Kokkos::fmin(Kokkos::fmax(v, 0.0), 1.0);
  const double w = plicAlphaReduced(Kokkos::fmin(vc, 1.0 - vc), n1, n2, n3);
  const double a1 = (vc <= 0.5) ? w : 1.0 - w;  // undo the V -> 1-V fold
 
  double alpha = s * a1;  // undo the L1 normalization
  if (mx < 0.0)
    alpha -= ax;  // undo the octant mirroring
  if (my < 0.0)
    alpha -= ay;
  if (mz < 0.0)
    alpha -= az;
  return alpha;
}
 
// ---------------------------------------------------------------------------------------------
// (5) slab truncation — the geometric face flux
// ---------------------------------------------------------------------------------------------
 
KOKKOS_INLINE_FUNCTION double plicSlabVolume(double mx, double my, double mz, double alpha, int dir,
                                             double a, double b) {
  const double len = b - a;
  if (!(len > 0.0))
    return 0.0;
  double m[3] = {mx, my, mz};
  const double md = m[dir];
  m[dir] = md * len;
  return len * plicVolume(m[0], m[1], m[2], alpha - md * a);
}
 
KOKKOS_INLINE_FUNCTION double plicBoxVolume(double mx, double my, double mz, double alpha,
                                            double a0, double b0, double a1, double b1, double a2,
                                            double b2) {
  const double len[3] = {b0 - a0, b1 - a1, b2 - a2};
  if (!(len[0] > 0.0) || !(len[1] > 0.0) || !(len[2] > 0.0))
    return 0.0;
  const double m[3] = {mx, my, mz};
  const double off = alpha - (m[0] * a0 + m[1] * a1 + m[2] * a2);
  return len[0] * len[1] * len[2] * plicVolume(m[0] * len[0], m[1] * len[1], m[2] * len[2], off);
}
 
KOKKOS_INLINE_FUNCTION double faceFluxVolume(double mx, double my, double mz, double alpha, int dir,
                                             double f) {
  return plicSlabVolume(mx, my, mz, alpha, dir, 0.0, f);
}
 
// ---------------------------------------------------------------------------------------------
// (3)/(4) interface normals from a 3x3x3 colour stencil
// ---------------------------------------------------------------------------------------------
 
KOKKOS_INLINE_FUNCTION int plicSt(int i, int j, int k) {
  return i + 3 * j + 9 * k;
}
 
#define PECLET_VOF_C(dx, dy, dz) c[plicSt(1 + (dx), 1 + (dy), 1 + (dz))]
 
KOKKOS_INLINE_FUNCTION void youngsNormal(const double c[27], double m[3]) {
  double m1, m2;
 
  m1 = PECLET_VOF_C(-1, -1, -1) + PECLET_VOF_C(-1, 1, -1) + PECLET_VOF_C(-1, -1, 1) +
       PECLET_VOF_C(-1, 1, 1) +
       2.0 * (PECLET_VOF_C(-1, -1, 0) + PECLET_VOF_C(-1, 1, 0) + PECLET_VOF_C(-1, 0, -1) +
              PECLET_VOF_C(-1, 0, 1)) +
       4.0 * PECLET_VOF_C(-1, 0, 0);
  m2 = PECLET_VOF_C(1, -1, -1) + PECLET_VOF_C(1, 1, -1) + PECLET_VOF_C(1, -1, 1) +
       PECLET_VOF_C(1, 1, 1) +
       2.0 * (PECLET_VOF_C(1, -1, 0) + PECLET_VOF_C(1, 1, 0) + PECLET_VOF_C(1, 0, -1) +
              PECLET_VOF_C(1, 0, 1)) +
       4.0 * PECLET_VOF_C(1, 0, 0);
  m[0] = m1 - m2;
 
  m1 = PECLET_VOF_C(-1, -1, -1) + PECLET_VOF_C(-1, -1, 1) + PECLET_VOF_C(1, -1, -1) +
       PECLET_VOF_C(1, -1, 1) +
       2.0 * (PECLET_VOF_C(-1, -1, 0) + PECLET_VOF_C(1, -1, 0) + PECLET_VOF_C(0, -1, -1) +
              PECLET_VOF_C(0, -1, 1)) +
       4.0 * PECLET_VOF_C(0, -1, 0);
  m2 = PECLET_VOF_C(-1, 1, -1) + PECLET_VOF_C(-1, 1, 1) + PECLET_VOF_C(1, 1, -1) +
       PECLET_VOF_C(1, 1, 1) +
       2.0 * (PECLET_VOF_C(-1, 1, 0) + PECLET_VOF_C(1, 1, 0) + PECLET_VOF_C(0, 1, -1) +
              PECLET_VOF_C(0, 1, 1)) +
       4.0 * PECLET_VOF_C(0, 1, 0);
  m[1] = m1 - m2;
 
  m1 = PECLET_VOF_C(-1, -1, -1) + PECLET_VOF_C(-1, 1, -1) + PECLET_VOF_C(1, -1, -1) +
       PECLET_VOF_C(1, 1, -1) +
       2.0 * (PECLET_VOF_C(-1, 0, -1) + PECLET_VOF_C(1, 0, -1) + PECLET_VOF_C(0, -1, -1) +
              PECLET_VOF_C(0, 1, -1)) +
       4.0 * PECLET_VOF_C(0, 0, -1);
  m2 = PECLET_VOF_C(-1, -1, 1) + PECLET_VOF_C(-1, 1, 1) + PECLET_VOF_C(1, -1, 1) +
       PECLET_VOF_C(1, 1, 1) +
       2.0 * (PECLET_VOF_C(-1, 0, 1) + PECLET_VOF_C(1, 0, 1) + PECLET_VOF_C(0, -1, 1) +
              PECLET_VOF_C(0, 1, 1)) +
       4.0 * PECLET_VOF_C(0, 0, 1);
  m[2] = m1 - m2;
 
  if (plicNormalizeL1(m) < 1e-30) {
    m[0] = 1.0;
    m[1] = 0.0;
    m[2] = 0.0;
  }
}
 
KOKKOS_INLINE_FUNCTION void mycNormal(const double c[27], double m[3]) {
  double m1, m2, mm[4][3], t0, t1, t2;
  int cn;
 
  // plane as sgn(mx) X = m01 Y + m02 Z + alpha
  m1 = PECLET_VOF_C(-1, 0, -1) + PECLET_VOF_C(-1, 0, 1) + PECLET_VOF_C(-1, -1, 0) +
       PECLET_VOF_C(-1, 1, 0) + PECLET_VOF_C(-1, 0, 0);
  m2 = PECLET_VOF_C(1, 0, -1) + PECLET_VOF_C(1, 0, 1) + PECLET_VOF_C(1, -1, 0) +
       PECLET_VOF_C(1, 1, 0) + PECLET_VOF_C(1, 0, 0);
  mm[0][0] = m1 > m2 ? 1.0 : -1.0;
 
  m1 = PECLET_VOF_C(-1, -1, 0) + PECLET_VOF_C(1, -1, 0) + PECLET_VOF_C(0, -1, 0);
  m2 = PECLET_VOF_C(-1, 1, 0) + PECLET_VOF_C(1, 1, 0) + PECLET_VOF_C(0, 1, 0);
  mm[0][1] = 0.5 * (m1 - m2);
 
  m1 = PECLET_VOF_C(-1, 0, -1) + PECLET_VOF_C(1, 0, -1) + PECLET_VOF_C(0, 0, -1);
  m2 = PECLET_VOF_C(-1, 0, 1) + PECLET_VOF_C(1, 0, 1) + PECLET_VOF_C(0, 0, 1);
  mm[0][2] = 0.5 * (m1 - m2);
 
  // plane as sgn(my) Y = m10 X + m12 Z + alpha
  m1 = PECLET_VOF_C(-1, -1, 0) + PECLET_VOF_C(-1, 1, 0) + PECLET_VOF_C(-1, 0, 0);
  m2 = PECLET_VOF_C(1, -1, 0) + PECLET_VOF_C(1, 1, 0) + PECLET_VOF_C(1, 0, 0);
  mm[1][0] = 0.5 * (m1 - m2);
 
  m1 = PECLET_VOF_C(0, -1, -1) + PECLET_VOF_C(0, -1, 1) + PECLET_VOF_C(1, -1, 0) +
       PECLET_VOF_C(-1, -1, 0) + PECLET_VOF_C(0, -1, 0);
  m2 = PECLET_VOF_C(0, 1, -1) + PECLET_VOF_C(0, 1, 1) + PECLET_VOF_C(1, 1, 0) +
       PECLET_VOF_C(-1, 1, 0) + PECLET_VOF_C(0, 1, 0);
  mm[1][1] = m1 > m2 ? 1.0 : -1.0;
 
  m1 = PECLET_VOF_C(0, -1, -1) + PECLET_VOF_C(0, 0, -1) + PECLET_VOF_C(0, 1, -1);
  m2 = PECLET_VOF_C(0, -1, 1) + PECLET_VOF_C(0, 0, 1) + PECLET_VOF_C(0, 1, 1);
  mm[1][2] = 0.5 * (m1 - m2);
 
  // plane as sgn(mz) Z = m20 X + m21 Y + alpha
  m1 = PECLET_VOF_C(-1, 0, -1) + PECLET_VOF_C(-1, 0, 1) + PECLET_VOF_C(-1, 0, 0);
  m2 = PECLET_VOF_C(1, 0, -1) + PECLET_VOF_C(1, 0, 1) + PECLET_VOF_C(1, 0, 0);
  mm[2][0] = 0.5 * (m1 - m2);
 
  m1 = PECLET_VOF_C(0, -1, -1) + PECLET_VOF_C(0, -1, 1) + PECLET_VOF_C(0, -1, 0);
  m2 = PECLET_VOF_C(0, 1, -1) + PECLET_VOF_C(0, 1, 1) + PECLET_VOF_C(0, 1, 0);
  mm[2][1] = 0.5 * (m1 - m2);
 
  m1 = PECLET_VOF_C(-1, 0, -1) + PECLET_VOF_C(1, 0, -1) + PECLET_VOF_C(0, -1, -1) +
       PECLET_VOF_C(0, 1, -1) + PECLET_VOF_C(0, 0, -1);
  m2 = PECLET_VOF_C(-1, 0, 1) + PECLET_VOF_C(1, 0, 1) + PECLET_VOF_C(0, -1, 1) +
       PECLET_VOF_C(0, 1, 1) + PECLET_VOF_C(0, 0, 1);
  mm[2][2] = m1 > m2 ? 1.0 : -1.0;
 
  plicNormalizeL1(mm[0]);
  plicNormalizeL1(mm[1]);
  plicNormalizeL1(mm[2]);
 
  // pick the centred candidate with the largest dominant component
  t0 = Kokkos::fabs(mm[0][0]);
  t1 = Kokkos::fabs(mm[1][1]);
  t2 = Kokkos::fabs(mm[2][2]);
  cn = 0;
  if (t1 > t0) {
    t0 = t1;
    cn = 1;
  }
  if (t2 > t0)
    cn = 2;
 
  youngsNormal(c, mm[3]);  // already L1-normalized; (1,0,0) on a degenerate stencil
 
  t0 = Kokkos::fabs(mm[3][0]);
  t1 = Kokkos::fabs(mm[3][1]);
  t2 = Kokkos::fabs(mm[3][2]);
  if (t1 > t0)
    t0 = t1;
  if (t2 > t0)
    t0 = t2;
  if (Kokkos::fabs(mm[cn][cn]) > t0)
    cn = 3;
 
  m[0] = mm[cn][0];
  m[1] = mm[cn][1];
  m[2] = mm[cn][2];
}
 
#undef PECLET_VOF_C
 
// ---------------------------------------------------------------------------------------------
// (6) test helpers — exact per-cell fractions for manufactured interfaces
// ---------------------------------------------------------------------------------------------
 
KOKKOS_INLINE_FUNCTION double planeCellFraction(double mx, double my, double mz, double alphaGlobal,
                                                double x0, double y0, double z0, double h) {
  const double aLoc = alphaGlobal - (mx * x0 + my * y0 + mz * z0);
  return plicVolume(mx * h, my * h, mz * h, aLoc);
}
 
KOKKOS_INLINE_FUNCTION double planeCellFractionAniso(double mx, double my, double mz,
                                                     double alphaGlobal, double x0, double y0,
                                                     double z0, const VofMetric& g) {
  const double aLoc = alphaGlobal - (mx * x0 + my * y0 + mz * z0);
  return plicVolume(mx * g.h[0], my * g.h[1], mz * g.h[2], aLoc);
}
 
KOKKOS_INLINE_FUNCTION double sphereCellFraction(double cx, double cy, double cz, double r,
                                                 double x0, double y0, double z0, double h,
                                                 int levels) {
  if (levels < 0)
    levels = 0;
  if (levels > 8)
    levels = 8;
  struct Box {
    double x, y, z, s;
    int lev;
  };
  Box stack[64];  // DFS over an octree of depth L keeps at most 7L+1 boxes live
  int top = 0;
  stack[top++] = Box{x0, y0, z0, h, 0};
  double vol = 0.0;
  const double kHalfDiag = 0.86602540378443865;  // sqrt(3)/2
 
  while (top > 0) {
    const Box b = stack[--top];
    const double half = 0.5 * b.s;
    const double dx = b.x + half - cx, dy = b.y + half - cy, dz = b.z + half - cz;
    const double d = Kokkos::sqrt(dx * dx + dy * dy + dz * dz);
    const double hd = kHalfDiag * b.s;
    if (r - d >= hd) {  // wholly inside the sphere
      vol += b.s * b.s * b.s;
      continue;
    }
    if (d - r >= hd)
      continue;  // wholly outside
    if (b.lev >= levels || !(d > 0.0)) {
      if (d > 0.0) {
        const double inv = 1.0 / d;
        const double nx = dx * inv, ny = dy * inv, nz = dz * inv;
        // fluid where (X - centre).n < r, i.e. n.X < r + n.centre
        vol += b.s * b.s * b.s *
               planeCellFraction(nx, ny, nz, r + nx * cx + ny * cy + nz * cz, b.x, b.y, b.z, b.s);
      }
      continue;
    }
    for (int q = 0; q < 8; ++q)
      stack[top++] = Box{b.x + ((q & 1) ? half : 0.0), b.y + ((q & 2) ? half : 0.0),
                         b.z + ((q & 4) ? half : 0.0), half, b.lev + 1};
  }
  return vol / (h * h * h);
}
 
}  // namespace peclet::core::vof
 
#endif  // PECLET_CORE_VOF_PLIC_HPP