marching_cubes.hpp

Go to the documentation of this file.

 
#ifndef PECLET_FLOW_VOF_MARCHING_CUBES_HPP
#define PECLET_FLOW_VOF_MARCHING_CUBES_HPP
 
#include <Kokkos_Core.hpp>
#include <Kokkos_MathematicalFunctions.hpp>
 
namespace peclet::flow::vof {
 
enum McSource : int {
  kMcSrcColour = 0,  
  kMcSrcPlic = 1,    
};
 
enum McDeposit : int {
  kMcDepositCentroid = 0,  
  kMcDepositSplit = 1,     
};
 
struct McVertex {
  double psi = 0.0;
  double d = 0.0;
  double n[3] = {0.0, 0.0, 0.0};
  bool has = false;
};
 
KOKKOS_INLINE_FUNCTION void mcCornerPos(int k, double p[3]) {
  p[0] = static_cast<double>(k & 1);
  p[1] = static_cast<double>((k >> 1) & 1);
  p[2] = static_cast<double>((k >> 2) & 1);
}
 
KOKKOS_INLINE_FUNCTION double mcEdgeT(const McVertex& a, const McVertex& b, const double dv[3],
                                      int src) {
  if (src == kMcSrcPlic) {
    double acc = 0.0;
    int cnt = 0;
    if (a.has) {
      const double g = a.n[0] * dv[0] + a.n[1] * dv[1] + a.n[2] * dv[2];
      if (Kokkos::fabs(g) > 1e-12) {
        const double t = -a.d / g;
        if (t > -0.5 && t < 1.5) {
          acc += t;
          ++cnt;
        }
      }
    }
    if (b.has) {
      const double g = b.n[0] * dv[0] + b.n[1] * dv[1] + b.n[2] * dv[2];
      if (Kokkos::fabs(g) > 1e-12) {
        const double t = 1.0 - b.d / g;
        if (t > -0.5 && t < 1.5) {
          acc += t;
          ++cnt;
        }
      }
    }
    if (cnt > 0) {
      double t = acc / static_cast<double>(cnt);
      if (a.has && b.has) {
        // Phi(s) = (1-s) phi_a(s) + s phi_b(s): the linear BLEND of the two cells' own signed
        // distance functions, which is exact on a plane and weights each cell's plane by how near
        // the crossing is to it (a quadratic in s).
        const double ga = a.n[0] * dv[0] + a.n[1] * dv[1] + a.n[2] * dv[2];
        const double gb = b.n[0] * dv[0] + b.n[1] * dv[1] + b.n[2] * dv[2];
        const double A = gb - ga, B = (ga - gb) + (b.d - a.d), C0 = a.d;
        double s2 = -1.0;
        if (Kokkos::fabs(A) < 1e-12 * (Kokkos::fabs(B) + 1e-30)) {
          if (Kokkos::fabs(B) > 0.0)
            s2 = -C0 / B;
        } else {
          const double disc = B * B - 4.0 * A * C0;
          if (disc >= 0.0) {
            const double sq = Kokkos::sqrt(disc);
            const double r1 = (-B + sq) / (2.0 * A), r2 = (-B - sq) / (2.0 * A);
            const bool o1 = (r1 >= 0.0 && r1 <= 1.0), o2 = (r2 >= 0.0 && r2 <= 1.0);
            if (o1 && !o2)
              s2 = r1;
            else if (o2 && !o1)
              s2 = r2;
            else if (o1 && o2)
              s2 = (Kokkos::fabs(r1 - t) < Kokkos::fabs(r2 - t)) ? r1 : r2;
          }
        }
        if (s2 >= 0.0 && s2 <= 1.0)
          t = s2;
      }
      return t < 0.0 ? 0.0 : (t > 1.0 ? 1.0 : t);
    }
  }
  const double den = a.psi - b.psi;
  if (!(Kokkos::fabs(den) > 0.0))
    return 0.5;
  const double t = a.psi / den;
  return t < 0.0 ? 0.0 : (t > 1.0 ? 1.0 : t);
}
 
KOKKOS_INLINE_FUNCTION double mcTwiceArea(const double p0[3], const double p1[3],
                                          const double p2[3], const VofMetric& g) {
  // Phase 3: the sheet is CUT in index space (every vertex sits on a lattice edge, which is what
  // makes it watertight) and MEASURED in physical space — the two edge vectors are mapped through
  // H before the cross product. `g.h = {1,1,1}` multiplies each by 1.0: today's numbers.
  const double ux = (p1[0] - p0[0]) * g.h[0], uy = (p1[1] - p0[1]) * g.h[1],
               uz = (p1[2] - p0[2]) * g.h[2];
  const double vx = (p2[0] - p0[0]) * g.h[0], vy = (p2[1] - p0[1]) * g.h[1],
               vz = (p2[2] - p0[2]) * g.h[2];
  const double cx = uy * vz - uz * vy, cy = uz * vx - ux * vz, cz = ux * vy - uy * vx;
  return Kokkos::sqrt(cx * cx + cy * cy + cz * cz);
}
 
KOKKOS_INLINE_FUNCTION double mcPolygonArea(const double p[8][3], int n, const VofMetric& g) {
  if (n < 3)
    return 0.0;
  // As `mcTwiceArea`: clipped in index space, measured in physical space. The shoelace vector of
  // the mapped polygon is `det(H) H^-1` applied to the index one (the cofactor map), which for a
  // diagonal H is the per-component product below — and exactly the identity at g.h = {1,1,1}.
  const double c0 = g.h[1] * g.h[2], c1 = g.h[2] * g.h[0], c2 = g.h[0] * g.h[1];
  double sx = 0.0, sy = 0.0, sz = 0.0;
  for (int i = 0; i < n; ++i) {
    const int j = (i + 1 == n) ? 0 : i + 1;
    sx += p[i][1] * p[j][2] - p[i][2] * p[j][1];
    sy += p[i][2] * p[j][0] - p[i][0] * p[j][2];
    sz += p[i][0] * p[j][1] - p[i][1] * p[j][0];
  }
  sx *= c0;
  sy *= c1;
  sz *= c2;
  return 0.5 * Kokkos::sqrt(sx * sx + sy * sy + sz * sz);
}
 
KOKKOS_INLINE_FUNCTION int mcClipHalf(double p[8][3], int n, int axis, double c, bool keepAbove) {
  double q[8][3];
  int m = 0;
  for (int i = 0; i < n; ++i) {
    const int j = (i + 1 == n) ? 0 : i + 1;
    const double si = keepAbove ? (p[i][axis] - c) : (c - p[i][axis]);
    const double sj = keepAbove ? (p[j][axis] - c) : (c - p[j][axis]);
    const bool ini = si >= 0.0, inj = sj >= 0.0;
    if (ini && m < 8) {
      q[m][0] = p[i][0];
      q[m][1] = p[i][1];
      q[m][2] = p[i][2];
      ++m;
    }
    if (ini != inj && m < 8) {
      const double den = si - sj;
      const double t = (Kokkos::fabs(den) > 0.0) ? si / den : 0.0;
      q[m][0] = p[i][0] + t * (p[j][0] - p[i][0]);
      q[m][1] = p[i][1] + t * (p[j][1] - p[i][1]);
      q[m][2] = p[i][2] + t * (p[j][2] - p[i][2]);
      ++m;
    }
  }
  for (int i = 0; i < m; ++i) {
    p[i][0] = q[i][0];
    p[i][1] = q[i][1];
    p[i][2] = q[i][2];
  }
  return m;
}
 
KOKKOS_INLINE_FUNCTION double mcClipToCorner(const double p0[3], const double p1[3],
                                             const double p2[3], int k, const VofMetric& g) {
  double poly[8][3];
  for (int d = 0; d < 3; ++d) {
    poly[0][d] = p0[d];
    poly[1][d] = p1[d];
    poly[2][d] = p2[d];
  }
  int n = 3;
  for (int d = 0; d < 3 && n >= 3; ++d)
    n = mcClipHalf(poly, n, d, 0.5, ((k >> d) & 1) != 0);
  return mcPolygonArea(poly, n, g);
}
 
KOKKOS_INLINE_FUNCTION int mcNearestMasked(const double p[3], int mask) {
  int best = -1;
  double bd = 0.0;
  for (int k = 0; k < 8; ++k) {
    if (!((mask >> k) & 1))
      continue;
    double q[3];
    mcCornerPos(k, q);
    const double dx = p[0] - q[0], dy = p[1] - q[1], dz = p[2] - q[2];
    const double d2 = dx * dx + dy * dy + dz * dz;
    if (best < 0 || d2 < bd) {
      best = k;
      bd = d2;
    }
  }
  return best;
}
 
KOKKOS_INLINE_FUNCTION int mcRetarget(int k, int mask) {
  if (mask == 0 || ((mask >> k) & 1))
    return k;
  double q[3];
  mcCornerPos(k, q);
  const int t = mcNearestMasked(q, mask);
  return t < 0 ? k : t;
}
 
KOKKOS_INLINE_FUNCTION double mcTriangleToCorner(const double p0[3], const double p1[3],
                                                 const double p2[3], int lc, int deposit, int mask,
                                                 const VofMetric& g) {
  if (deposit == kMcDepositCentroid) {
    double ctr[3];
    for (int d = 0; d < 3; ++d)
      ctr[d] = (p0[d] + p1[d] + p2[d]) * (1.0 / 3.0);
    int corner = 0;
    for (int d = 0; d < 3; ++d)
      if (ctr[d] >= 0.5)
        corner |= 1 << d;
    if (mask != 0 && !((mask >> corner) & 1)) {
      const int t = mcNearestMasked(ctr, mask);
      if (t >= 0)
        corner = t;
    }
    return (corner == lc) ? 0.5 * mcTwiceArea(p0, p1, p2, g) : 0.0;
  }
  double acc = 0.0;
  if (mask == 0 || ((mask >> lc) & 1))
    acc += mcClipToCorner(p0, p1, p2, lc, g);
  if (mask != 0 && mask != 255)
    for (int k = 0; k < 8; ++k)
      if (!((mask >> k) & 1) && mcRetarget(k, mask) == lc)
        acc += mcClipToCorner(p0, p1, p2, k, g);
  return acc;
}
 
KOKKOS_INLINE_FUNCTION double mcCubeCornerArea(const McVertex v[8], int lc, int src, int deposit,
                                               const VofMetric& g) {
  // Kuhn's decomposition along the body diagonal 0 -> 7. Corner index = x + 2y + 4z.
  const int TET[6][4] = {{0, 1, 3, 7}, {0, 1, 5, 7}, {0, 2, 3, 7},
                         {0, 2, 6, 7}, {0, 4, 5, 7}, {0, 4, 6, 7}};
  int mask = 0;
  for (int k = 0; k < 8; ++k)
    if (v[k].has)
      mask |= 1 << k;
  // No mixed cell in the cube = no reconstructible interface, and nothing this rung's consumer
  // could carry an `mdot` on. This is the same statement the rest of the rung makes (a grid-
  // aligned sharp interface has no interfacial cell and therefore no PLIC area either), and it is
  // what keeps a PERIODIC SEAM — the C = 1 | C = 0 jump a half-space wraps into — from being
  // reported as interface.
  if (mask == 0)
    return 0.0;
  double acc = 0.0;
  for (int t = 0; t < 6; ++t) {
    int code = 0;
    for (int j = 0; j < 4; ++j)
      if (v[TET[t][j]].psi < 0.0)  // liquid
        code |= 1 << j;
    if (code == 0 || code == 15)
      continue;
    int nneg = 0;
    for (int j = 0; j < 4; ++j)
      if (code & (1 << j))
        ++nneg;
    double tri[4][3];  // up to 4 crossing points (the 2-vs-2 quad)
    int ntri = 0;
    if (nneg == 1 || nneg == 3) {
      // the lone corner, and the three edges leaving it
      const bool lonesign = (nneg == 1);
      int lone = -1;
      for (int j = 0; j < 4; ++j)
        if (((code & (1 << j)) != 0) == lonesign) {
          lone = j;
          break;
        }
      for (int j = 0; j < 4; ++j) {
        if (j == lone)
          continue;
        const int ca = TET[t][lone], cb = TET[t][j];
        double pa[3], pb[3];
        mcCornerPos(ca, pa);
        mcCornerPos(cb, pb);
        const double dv[3] = {pb[0] - pa[0], pb[1] - pa[1], pb[2] - pa[2]};
        const double s = mcEdgeT(v[ca], v[cb], dv, src);
        for (int d = 0; d < 3; ++d)
          tri[ntri][d] = pa[d] + s * dv[d];
        ++ntri;
      }
      acc += mcTriangleToCorner(tri[0], tri[1], tri[2], lc, deposit, mask, g);
    } else {
      // 2 vs 2: the quad through the four edges that join the two pairs, in cyclic order
      int neg[2], pos[2], nn = 0, np = 0;
      for (int j = 0; j < 4; ++j)
        ((code & (1 << j)) ? neg[nn++] : pos[np++]) = j;
      const int pair[4][2] = {
          {neg[0], pos[0]}, {neg[0], pos[1]}, {neg[1], pos[1]}, {neg[1], pos[0]}};
      for (int q = 0; q < 4; ++q) {
        const int ca = TET[t][pair[q][0]], cb = TET[t][pair[q][1]];
        double pa[3], pb[3];
        mcCornerPos(ca, pa);
        mcCornerPos(cb, pb);
        const double dv[3] = {pb[0] - pa[0], pb[1] - pa[1], pb[2] - pa[2]};
        const double s = mcEdgeT(v[ca], v[cb], dv, src);
        for (int d = 0; d < 3; ++d)
          tri[q][d] = pa[d] + s * dv[d];
      }
      acc += mcTriangleToCorner(tri[0], tri[1], tri[2], lc, deposit, mask, g);
      acc += mcTriangleToCorner(tri[0], tri[2], tri[3], lc, deposit, mask, g);
    }
  }
  return acc;
}
 
KOKKOS_INLINE_FUNCTION double mcCubeCornerArea(const McVertex v[8], int lc, int src, int deposit) {
  return mcCubeCornerArea(v, lc, src, deposit, VofMetric{});
}
 
}  // namespace peclet::flow::vof
 
#endif  // PECLET_FLOW_VOF_MARCHING_CUBES_HPP