sdf.hpp

Go to the documentation of this file.

 
#ifndef PECLET_VORO_SDF_HPP
#define PECLET_VORO_SDF_HPP
 
#include <Kokkos_Core.hpp>
#include <type_traits>
 
#include "peclet/core/common/view.hpp"
#include "peclet/voro/convex_cell.hpp"
 
namespace peclet::voro {
 
constexpr int kBoundaryFacet = -2;
 
struct NoSdf {};
 
template <class Real>
struct SdfSphere {
  Real cx = 0, cy = 0, cz = 0, radius = 1;
  KOKKOS_INLINE_FUNCTION Real eval(Real x, Real y, Real z) const {
    Real dx = x - cx, dy = y - cy, dz = z - cz;
    return Kokkos::sqrt(dx * dx + dy * dy + dz * dz) - radius;
  }
  KOKKOS_INLINE_FUNCTION Real gradH() const { return Real(1e-4); }
};
 
template <class Real>
struct SdfBox {
  Real cx = 0, cy = 0, cz = 0, hx = Real(0.5), hy = Real(0.5), hz = Real(0.5);
  KOKKOS_INLINE_FUNCTION Real eval(Real x, Real y, Real z) const {
    Real qx = Kokkos::fabs(x - cx) - hx, qy = Kokkos::fabs(y - cy) - hy,
         qz = Kokkos::fabs(z - cz) - hz;
    Real ox = qx > 0 ? qx : 0, oy = qy > 0 ? qy : 0, oz = qz > 0 ? qz : 0;
    Real outside = Kokkos::sqrt(ox * ox + oy * oy + oz * oz);
    Real mx = qx > qy ? qx : qy;
    mx = mx > qz ? mx : qz;
    Real inside = mx < 0 ? mx : 0;
    return outside + inside;
  }
  KOKKOS_INLINE_FUNCTION Real gradH() const { return Real(1e-4); }
};
 
template <class Real>
struct SdfHollowCylinder {
  Real cx = 0, cy = 0, cz = 0, rOuter = 1, rInner = Real(0.5), height = 1;
  int axis = 2;
  KOKKOS_INLINE_FUNCTION Real eval(Real x, Real y, Real z) const {
    Real d[3] = {x - cx, y - cy, z - cz};
    int a0 = (axis + 1) % 3, a1 = (axis + 2) % 3;
    Real r = Kokkos::sqrt(d[a0] * d[a0] + d[a1] * d[a1]);
    Real zz = d[axis];
    Real m = r - rOuter;
    Real b = rInner - r;
    if (b > m)
      m = b;
    Real c = Kokkos::fabs(zz) - height * Real(0.5);
    if (c > m)
      m = c;
    return m;
  }
  KOKKOS_INLINE_FUNCTION Real gradH() const { return Real(1e-4); }
};
 
template <class Real>
struct SdfGrid {
  Kokkos::View<const float*, peclet::core::MemSpace> values;  // i + nx*(j + ny*k)
  int nx = 0, ny = 0, nz = 0;
  Real ox = 0, oy = 0, oz = 0, sx = 1, sy = 1, sz = 1;
  KOKKOS_INLINE_FUNCTION Real at(int i, int j, int k) const {
    return static_cast<Real>(values(i + nx * (j + ny * k)));
  }
  KOKKOS_INLINE_FUNCTION Real eval(Real x, Real y, Real z) const {
    Real p[3] = {x, y, z};
    Real o[3] = {ox, oy, oz}, s[3] = {sx, sy, sz};
    int dims[3] = {nx, ny, nz};
    int i0[3], i1[3];
    Real f[3];
    for (int d = 0; d < 3; ++d) {
      Real c = (p[d] - o[d]) / s[d];
      Real hi = static_cast<Real>(dims[d] - 1);
      c = c < 0 ? Real(0) : (c > hi ? hi : c);
      i0[d] = static_cast<int>(Kokkos::floor(c));
      i1[d] = i0[d] + 1 < dims[d] ? i0[d] + 1 : dims[d] - 1;
      f[d] = c - static_cast<Real>(i0[d]);
    }
    Real c000 = at(i0[0], i0[1], i0[2]), c100 = at(i1[0], i0[1], i0[2]);
    Real c010 = at(i0[0], i1[1], i0[2]), c110 = at(i1[0], i1[1], i0[2]);
    Real c001 = at(i0[0], i0[1], i1[2]), c101 = at(i1[0], i0[1], i1[2]);
    Real c011 = at(i0[0], i1[1], i1[2]), c111 = at(i1[0], i1[1], i1[2]);
    Real c00 = c000 * (1 - f[0]) + c100 * f[0];
    Real c10 = c010 * (1 - f[0]) + c110 * f[0];
    Real c01 = c001 * (1 - f[0]) + c101 * f[0];
    Real c11 = c011 * (1 - f[0]) + c111 * f[0];
    Real c0 = c00 * (1 - f[1]) + c10 * f[1];
    Real c1 = c01 * (1 - f[1]) + c11 * f[1];
    return c0 * (1 - f[2]) + c1 * f[2];
  }
  KOKKOS_INLINE_FUNCTION Real gradH() const {
    Real m = sx < sy ? sx : sy;
    m = m < sz ? m : sz;
    return Real(0.25) * m;
  }
};
 
template <class Real>
struct SdfSpheres {
  Kokkos::View<const Real*, peclet::core::MemSpace> cen;  // 3*M
  Kokkos::View<const Real*, peclet::core::MemSpace> rad;  // M
  int n = 0;
  Real L = 0;
  KOKKOS_INLINE_FUNCTION Real eval(Real x, Real y, Real z) const {
    Real m = Real(1e30);
    for (int i = 0; i < n; ++i) {
      Real dx = x - cen(3 * i), dy = y - cen(3 * i + 1), dz = z - cen(3 * i + 2);
      if (L > Real(0)) {
        dx -= L * Kokkos::round(dx / L);
        dy -= L * Kokkos::round(dy / L);
        dz -= L * Kokkos::round(dz / L);
      }
      const Real d = Kokkos::sqrt(dx * dx + dy * dy + dz * dz) - rad(i);
      if (d < m) m = d;
    }
    return m;
  }
  KOKKOS_INLINE_FUNCTION Real gradH() const { return Real(1e-4); }
};
 
template <class Real, class Sdf>
KOKKOS_INLINE_FUNCTION void sdfGradient(const Sdf& s, Real x, Real y, Real z, Real g[3]) {
  const Real h = s.gradH();
  g[0] = (s.eval(x + h, y, z) - s.eval(x - h, y, z)) / (2 * h);
  g[1] = (s.eval(x, y + h, z) - s.eval(x, y - h, z)) / (2 * h);
  g[2] = (s.eval(x, y, z + h) - s.eval(x, y, z - h)) / (2 * h);
}
 
template <class Real, class Sdf>
KOKKOS_INLINE_FUNCTION void sdfHessian(const Sdf& s, Real x, Real y, Real z, Real H[3][3]) {
  const Real h = s.gradH();
  Real gp[3], gm[3];
  sdfGradient<Real>(s, x + h, y, z, gp);
  sdfGradient<Real>(s, x - h, y, z, gm);
  for (int r = 0; r < 3; ++r) H[r][0] = (gp[r] - gm[r]) / (2 * h);
  sdfGradient<Real>(s, x, y + h, z, gp);
  sdfGradient<Real>(s, x, y - h, z, gm);
  for (int r = 0; r < 3; ++r) H[r][1] = (gp[r] - gm[r]) / (2 * h);
  sdfGradient<Real>(s, x, y, z + h, gp);
  sdfGradient<Real>(s, x, y, z - h, gm);
  for (int r = 0; r < 3; ++r) H[r][2] = (gp[r] - gm[r]) / (2 * h);
  for (int a = 0; a < 3; ++a)
    for (int b = a + 1; b < 3; ++b) {
      const Real m = Real(0.5) * (H[a][b] + H[b][a]);
      H[a][b] = H[b][a] = m;
    }
}
 
template <class Real, int MAXP, int MAXT, bool TrackAdj, class Sdf>
KOKKOS_INLINE_FUNCTION void addSdfWallForce(const ConvexCell<Real, MAXP, MAXT, TrackAdj>& c,
                                            const Real seed[3], const Sdf& sdf, const Real* gx,
                                            const Real* gy, const Real* gz, Real fSelf[3]) {
  if constexpr (std::is_same_v<Sdf, NoSdf>) {
    (void)c;
    (void)seed;
    (void)sdf;
    (void)gx;
    (void)gy;
    (void)gz;
    (void)fSelf;
    return;
  } else {
    // Aggregate the wall facets' geometry gradients (all move together under the seed-foot model).
    Real gw[3] = {Real(0), Real(0), Real(0)};
    bool any = false;
    for (int k = 0; k < c.np; ++k)
      if (c.pnbr[k] == kBoundaryFacet) {
        gw[0] += gx[k];
        gw[1] += gy[k];
        gw[2] += gz[k];
        any = true;
      }
    if (!any) return;
    const Real phi = sdf.eval(seed[0], seed[1], seed[2]);
    Real g[3];
    sdfGradient<Real>(sdf, seed[0], seed[1], seed[2], g);
    const Real gn = Kokkos::sqrt(g[0] * g[0] + g[1] * g[1] + g[2] * g[2]);
    if (gn <= Real(1e-12)) return;  // degenerate gradient (box crease/corner) — no wall force (guard)
    const Real u[3] = {g[0] / gn, g[1] / gn, g[2] / gn};
    Real H[3][3];
    sdfHessian<Real>(sdf, seed[0], seed[1], seed[2], H);
    const Real ug = u[0] * gw[0] + u[1] * gw[1] + u[2] * gw[2];
    const Real perp[3] = {gw[0] - ug * u[0], gw[1] - ug * u[1], gw[2] - ug * u[2]};
    const Real Hp[3] = {H[0][0] * perp[0] + H[0][1] * perp[1] + H[0][2] * perp[2],
                        H[1][0] * perp[0] + H[1][1] * perp[1] + H[1][2] * perp[2],
                        H[2][0] * perp[0] + H[2][1] * perp[1] + H[2][2] * perp[2]};
    const Real cH = phi / gn;
    for (int d = 0; d < 3; ++d) fSelf[d] += -gn * ug * u[d] - cH * Hp[d];
  }
}
 
template <class Real, int MAXP, int MAXT, bool TrackAdj, class Sdf>
KOKKOS_INLINE_FUNCTION bool clipCellAgainstSdf(ConvexCell<Real, MAXP, MAXT, TrackAdj>& c,
                                               const Real seed[3], const Sdf& sdf) {
  const Real tol = Real(1e-8);
  const int maxCuts = 24;
  const Real phiCenter = sdf.eval(seed[0], seed[1], seed[2]);
  if (phiCenter <= Real(0)) {  // seed inside solid -> no cell
    for (int t = 0; t < c.nt; ++t)
      c.alive[t] = false;
    return true;
  }
  // cell circumradius (dual vertices are seed-relative)
  const Real maxRsq = c.maxVertexRsq();
  const Real radius = Kokkos::sqrt(maxRsq > 0 ? maxRsq : Real(0));
  if (phiCenter > radius + tol)
    return false;  // cell fully in fluid
 
  bool seedPlaneApplied = false;
  for (int iter = 0; iter < maxCuts; ++iter) {
    Real probe[3] = {seed[0], seed[1], seed[2]};
    Real probePhi = phiCenter;
    if (seedPlaneApplied) {
      bool found = false;
      for (int t = 0; t < c.nt; ++t) {
        if (!c.alive[t])
          continue;
        Real x = seed[0] + c.vx[t], y = seed[1] + c.vy[t], z = seed[2] + c.vz[t];
        Real phi = sdf.eval(x, y, z);
        if (!found || phi < probePhi) {
          probe[0] = x;
          probe[1] = y;
          probe[2] = z;
          probePhi = phi;
          found = true;
        }
      }
      if (!found || probePhi >= -tol)
        break;
    }
    // closest point on sdf=0 + outward normal
    Real g[3];
    sdfGradient<Real>(sdf, probe[0], probe[1], probe[2], g);
    Real gsq = g[0] * g[0] + g[1] * g[1] + g[2] * g[2];
    if (gsq <= Real(0))
      break;
    Real phi = sdf.eval(probe[0], probe[1], probe[2]);
    Real invGsq = Real(1) / gsq, invG = Real(1) / Kokkos::sqrt(gsq);
    Real surf[3], normal[3];
    for (int k = 0; k < 3; ++k) {
      surf[k] = probe[k] - phi * g[k] * invGsq;
      normal[k] = g[k] * invG;
    }
    // orient the normal into the fluid (positive side)
    Real eps = Real(1e-3) * (radius + Real(1));
    if (eps < Real(1e-6))
      eps = Real(1e-6);
    if (sdf.eval(surf[0] + eps * normal[0], surf[1] + eps * normal[1], surf[2] + eps * normal[2]) <=
        Real(0))
      for (int k = 0; k < 3; ++k)
        normal[k] = -normal[k];
    // plane in seed-relative coords: dist(v) = v·pv - off; vertices in the solid (off the fluid
    // side) get dist > 0 and are removed.
    Real pv[3] = {-normal[0], -normal[1], -normal[2]};
    Real off =
        pv[0] * (surf[0] - seed[0]) + pv[1] * (surf[1] - seed[1]) + pv[2] * (surf[2] - seed[2]);
    c.clip(pv, off, kBoundaryFacet);
    seedPlaneApplied = true;
    if (c.empty())
      break;
  }
  return false;
}
 
}  // namespace peclet::voro
 
#endif  // PECLET_VORO_SDF_HPP