pore_cells.hpp

Go to the documentation of this file.

 
#ifndef PECLET_VORO_PORE_CELLS_HPP
#define PECLET_VORO_PORE_CELLS_HPP
 
#include <cstdint>
#include <Kokkos_Core.hpp>
#include <string>
 
#include "peclet/core/common/view.hpp"
#include "peclet/voro/convex_cell.hpp"
#include "peclet/voro/params.hpp"
#include "peclet/voro/sdf.hpp"
#include "peclet/voro/tess_grid.hpp"
#include "peclet/voro/tessellator.hpp"
 
namespace peclet::voro {
 
template <class Real>
struct PoreCellsResult {
  using MemSpace = peclet::core::MemSpace;
  Kokkos::View<Real*, MemSpace> points;  
  Kokkos::View<int64_t*, MemSpace> faces;
  Kokkos::View<int64_t*, MemSpace> faceOffset;  
  Kokkos::View<Real*, MemSpace> volume;         
  Kokkos::View<int*, MemSpace> boundary;        
  Kokkos::View<int*, MemSpace> seed;            
  long numOverflow = 0;                         
  long numIncomplete = 0;  
};
 
template <class Real>
struct PoreSectionResult {
  using MemSpace = peclet::core::MemSpace;
  Kokkos::View<Real*, MemSpace> verts;      
  Kokkos::View<int64_t*, MemSpace> offset;  
  Kokkos::View<Real*, MemSpace> volume;     
  Kokkos::View<int*, MemSpace> seed;        
  long numOverflow = 0;                     
  long numIncomplete = 0;                   
};
 
namespace detail {
 
template <class Cell>
KOKKOS_INLINE_FUNCTION bool poreCellCount(const Cell& c, int& nPts, int& nEntries, bool& wall) {
  if (c.empty() || c.overflow)
    return false;
  nPts = 0;
  for (int t = 0; t < c.nt; ++t)
    if (c.alive[t])
      ++nPts;
  if (nPts < 4)
    return false;
  nEntries = 1;
  wall = false;
  int nf = 0;
  int fidx[Cell::MAXFV];
  for (int k = 0; k < c.np; ++k) {
    const int m = c.faceOrderedIdx(k, fidx);
    if (m < 3)
      continue;
    ++nf;
    nEntries += 1 + m;
    if (c.pnbr[k] == kBoundaryFacet)
      wall = true;
  }
  return nf >= 4;
}
 
template <class Cell, class Real>
KOKKOS_INLINE_FUNCTION void poreCellFill(const Cell& c, const Real seed[3], int64_t ptBase,
                                         int64_t fBase, Real* pts, int64_t* faces) {
  int triToPt[Cell::kMaxT];
  int q = 0;
  for (int t = 0; t < c.nt; ++t) {
    triToPt[t] = -1;
    if (!c.alive[t])
      continue;
    triToPt[t] = q;
    const int64_t o = 3 * (ptBase + q);
    pts[o] = seed[0] + c.vx[t];
    pts[o + 1] = seed[1] + c.vy[t];
    pts[o + 2] = seed[2] + c.vz[t];
    ++q;
  }
  int fidx[Cell::MAXFV];
  int nf = 0;
  for (int k = 0; k < c.np; ++k)
    if (c.faceOrderedIdx(k, fidx) >= 3)
      ++nf;
  int64_t w = fBase;
  faces[w++] = nf;
  for (int k = 0; k < c.np; ++k) {
    const int m = c.faceOrderedIdx(k, fidx);
    if (m < 3)
      continue;
    faces[w++] = m;
    for (int i = 0; i < m; ++i)
      faces[w++] = ptBase + triToPt[fidx[i]];
  }
}
 
template <class Real, class Sdf, int MAXP, int MAXT>
struct PoreCellGather {
  using Builder = CellBuilder<Real, false, Sdf, false, MAXP, MAXT>;
  using Cell = typename Builder::Cell;
  Builder op;
 
  static PoreCellGather make(const TessGrid<Real>& grid, const Sdf& sdf) {
    PoreCellGather g{};
    Builder& op = g.op;
    op.binned = grid.binned;
    op.posSorted = grid.posSorted;
    op.wSorted = grid.wSorted;
    op.gidSorted = grid.gidSorted;
    op.cellStart = grid.cellStart;
    op.wlOff = grid.wlOff;
    op.wlRmin = grid.wlRmin;
    op.icx = grid.icx;
    op.icy = grid.icy;
    op.icz = grid.icz;
    op.Lx = grid.Lx;
    op.Ly = grid.Ly;
    op.Lz = grid.Lz;
    op.minCsz = grid.minCsz;
    op.wMaxAll = Real(0);
    op.dimx = grid.dimx;
    op.dimy = grid.dimy;
    op.dimz = grid.dimz;
    op.sw = grid.sw;
    op.nOff = grid.nOff;
    op.wlS = grid.wlS;
    op.useMorton = grid.useMorton;
    op.haveGid = grid.haveGid;
    op.sdf = sdf;
    return g;
  }
 
  KOKKOS_INLINE_FUNCTION int clipped(int pi, Cell& c, Real seedW[3]) const {
    Real wSelf;
    int ncRec, nnRec;
    op.gatherCell(pi, c, seedW[0], seedW[1], seedW[2], wSelf, ncRec, nnRec);
    if (c.overflow)
      return kOverflow;
    clipCellAgainstSdf<Real, MAXP, MAXT>(c, seedW, op.sdf);
    if (c.overflow)
      return kOverflow;
    if (c.empty())
      return kEmpty;
    // Completeness of the EXPORTED (clipped) cell: every seed that could cut it lies within twice
    // its reach, and the worklist examined everything within the walked window (grid.sw blocks,
    // the clamped window the grid was actually built with). The cold build judges the un-clipped
    // cell, which near a solid extends across the hole the clip removes — a spurious flag here.
    const Real covSq = Real(op.sw) * op.minCsz * Real(op.sw) * op.minCsz;
    return (covSq > Voronoi::template blockReachSq<Real>(c.maxVertexRsq(), wSelf, op.wMaxAll))
               ? kOk
               : kIncomplete;
  }
};
 
template <class Count, class Off>
inline void exclusiveScan(const Count& count, const Off& off, int N) {
  using Exec = peclet::core::ExecSpace;
  using T = typename Off::value_type;
  Kokkos::parallel_scan(
      "pore.scan", Kokkos::RangePolicy<Exec>(0, N + 1),
      KOKKOS_LAMBDA(const int i, T& acc, const bool fin) {
        const T v = (i < N) ? T(count(i)) : T(0);
        if (fin)
          off(i) = acc;
        acc += v;
      });
}
 
template <class Off>
inline typename Off::value_type lastOf(const Off& off, int N) {
  typename Off::value_type v;
  Kokkos::deep_copy(v, Kokkos::subview(off, N));
  return v;
}
 
template <class Status>
inline void statusCounts(const Status& status, int N, long& overflow, long& incomplete) {
  using Exec = peclet::core::ExecSpace;
  long ov = 0, inc = 0;
  Kokkos::parallel_reduce(
      "pore.overflow", Kokkos::RangePolicy<Exec>(0, N),
      KOKKOS_LAMBDA(const int i, long& a) { a += (status(i) & kOverflow) ? 1 : 0; }, ov);
  Kokkos::parallel_reduce(
      "pore.incomplete", Kokkos::RangePolicy<Exec>(0, N),
      KOKKOS_LAMBDA(const int i, long& a) {
        a += ((status(i) & kIncomplete) && !(status(i) & (kEmpty | kOverflow))) ? 1 : 0;
      },
      inc);
  overflow = ov;
  incomplete = inc;
}
 
}  // namespace detail
 
template <class Real, class Sdf, int MAXP = kMaxPlanes, int MAXT = kMaxTriangles>
PoreCellsResult<Real> buildPoreCells(const Kokkos::View<Real*, peclet::core::MemSpace>& posFlat,
                                     int N, const Real L[3], const Sdf& sdf,
                                     int sw = kSearchWindow) {
  using MemSpace = peclet::core::MemSpace;
  using Exec = peclet::core::ExecSpace;
  using Gather = detail::PoreCellGather<Real, Sdf, MAXP, MAXT>;
  using Cell = typename Gather::Cell;
  using Kokkos::view_alloc;
  using Kokkos::WithoutInitializing;
  Kokkos::View<Real*, MemSpace> noW;
  Kokkos::View<long*, MemSpace> noGid;
  auto grid = buildTessGrid<Real, false>(posFlat, noW, N, L, sw, N, noGid, nullptr);
  const Gather g = Gather::make(grid, sdf);
  auto binned = grid.binned;
 
  // Count pass (written at the ORIGINAL seed index, so the scan below is in seed order).
  Kokkos::View<int*, MemSpace> status("pore.status", N), emit("pore.emit", N), nPts("pore.nPts", N),
      nEnt("pore.nEnt", N), wall("pore.wall", N);
  Kokkos::parallel_for(
      "pore.count", Kokkos::RangePolicy<Exec>(0, N), KOKKOS_LAMBDA(const int pi) {
        const int i = binned(pi);
        Cell c;
        Real s[3];
        const int st = g.clipped(pi, c, s);
        status(i) = st;
        if (st & (kEmpty | kOverflow))
          return;
        int np = 0, ne = 0;
        bool w = false;
        if (!detail::poreCellCount(c, np, ne, w))
          return;
        emit(i) = 1;
        nPts(i) = np;
        nEnt(i) = ne;
        wall(i) = w ? 1 : 0;
      });
  Kokkos::View<int*, MemSpace> cellOff("pore.cellOff", N + 1);
  Kokkos::View<int64_t*, MemSpace> ptOff("pore.ptOff", N + 1), entOff("pore.entOff", N + 1);
  detail::exclusiveScan(emit, cellOff, N);
  detail::exclusiveScan(nPts, ptOff, N);
  detail::exclusiveScan(nEnt, entOff, N);
  const int nCells = detail::lastOf(cellOff, N);
  const int64_t nPoints = detail::lastOf(ptOff, N), nEntries = detail::lastOf(entOff, N);
 
  PoreCellsResult<Real> res;
  res.points = Kokkos::View<Real*, MemSpace>(
      view_alloc(std::string("pore.points"), WithoutInitializing), 3 * (size_t)nPoints);
  res.faces = Kokkos::View<int64_t*, MemSpace>(
      view_alloc(std::string("pore.faces"), WithoutInitializing), (size_t)nEntries);
  res.faceOffset = Kokkos::View<int64_t*, MemSpace>("pore.faceOffset", nCells + 1);
  res.volume = Kokkos::View<Real*, MemSpace>("pore.volume", nCells);
  res.boundary = Kokkos::View<int*, MemSpace>("pore.boundary", nCells);
  res.seed = Kokkos::View<int*, MemSpace>("pore.seed", nCells);
  auto pts = res.points;
  auto faces = res.faces;
  auto faceOffset = res.faceOffset;
  auto volume = res.volume;
  auto boundary = res.boundary;
  auto seed = res.seed;
  Kokkos::deep_copy(Kokkos::subview(faceOffset, nCells), nEntries);
  // Fill pass: the same gather again, written at the scanned offsets.
  Kokkos::parallel_for(
      "pore.fill", Kokkos::RangePolicy<Exec>(0, N), KOKKOS_LAMBDA(const int pi) {
        const int i = binned(pi);
        if (!emit(i))
          return;
        Cell c;
        Real s[3];
        g.clipped(pi, c, s);
        const int ci = cellOff(i);
        detail::poreCellFill(c, s, ptOff(i), entOff(i), pts.data(), faces.data());
        faceOffset(ci) = entOff(i);
        volume(ci) = c.volumePerVertex();
        boundary(ci) = wall(i);
        seed(ci) = i;
      });
  Kokkos::fence();
  detail::statusCounts(status, N, res.numOverflow, res.numIncomplete);
  return res;
}
 
template <class Real, class Sdf, int MAXP = kMaxPlanes, int MAXT = kMaxTriangles>
PoreSectionResult<Real> buildPoreSection(const Kokkos::View<Real*, peclet::core::MemSpace>& posFlat,
                                         int N, const Real L[3], const Sdf& sdf,
                                         const Real point[3], const Real normal[3],
                                         int sw = kSearchWindow) {
  using MemSpace = peclet::core::MemSpace;
  using Exec = peclet::core::ExecSpace;
  using Gather = detail::PoreCellGather<Real, Sdf, MAXP, MAXT>;
  using Cell = typename Gather::Cell;
  using Kokkos::view_alloc;
  using Kokkos::WithoutInitializing;
  Kokkos::View<Real*, MemSpace> noW;
  Kokkos::View<long*, MemSpace> noGid;
  auto grid = buildTessGrid<Real, false>(posFlat, noW, N, L, sw, N, noGid, nullptr);
  const Gather g = Gather::make(grid, sdf);
  auto binned = grid.binned;
  const Real px = point[0], py = point[1], pz = point[2];
  const Real ux = normal[0], uy = normal[1], uz = normal[2];
 
  Kokkos::View<int*, MemSpace> status("pore.status", N), nV("pore.nV", N), emit("pore.emit", N);
  Kokkos::parallel_for(
      "pore.section.count", Kokkos::RangePolicy<Exec>(0, N), KOKKOS_LAMBDA(const int pi) {
        const int i = binned(pi);
        Cell c;
        Real s[3];
        const int st = g.clipped(pi, c, s);
        status(i) = st;
        if (st & (kEmpty | kOverflow))
          return;
        const Real p0[3] = {px - s[0], py - s[1], pz - s[2]};  // the plane in the cell frame
        const Real u[3] = {ux, uy, uz};
        Real spx[Cell::MAXSV], spy[Cell::MAXSV], spz[Cell::MAXSV];
        const int m = c.sectionPolygon(p0, u, spx, spy, spz);
        if (m < 3)
          return;
        emit(i) = 1;
        nV(i) = m;
      });
  Kokkos::View<int*, MemSpace> polyOff("pore.polyOff", N + 1);
  Kokkos::View<int64_t*, MemSpace> vOff("pore.vOff", N + 1);
  detail::exclusiveScan(emit, polyOff, N);
  detail::exclusiveScan(nV, vOff, N);
  const int nPoly = detail::lastOf(polyOff, N);
  const int64_t nVerts = detail::lastOf(vOff, N);
 
  PoreSectionResult<Real> res;
  res.verts = Kokkos::View<Real*, MemSpace>(
      view_alloc(std::string("pore.verts"), WithoutInitializing), 3 * (size_t)nVerts);
  res.offset = Kokkos::View<int64_t*, MemSpace>("pore.offset", nPoly + 1);
  res.volume = Kokkos::View<Real*, MemSpace>("pore.volume", nPoly);
  res.seed = Kokkos::View<int*, MemSpace>("pore.seed", nPoly);
  auto verts = res.verts;
  auto offset = res.offset;
  auto volume = res.volume;
  auto seed = res.seed;
  Kokkos::deep_copy(Kokkos::subview(offset, nPoly), nVerts);
  Kokkos::parallel_for(
      "pore.section.fill", Kokkos::RangePolicy<Exec>(0, N), KOKKOS_LAMBDA(const int pi) {
        const int i = binned(pi);
        if (!emit(i))
          return;
        Cell c;
        Real s[3];
        g.clipped(pi, c, s);
        const Real p0[3] = {px - s[0], py - s[1], pz - s[2]};
        const Real u[3] = {ux, uy, uz};
        Real spx[Cell::MAXSV], spy[Cell::MAXSV], spz[Cell::MAXSV];
        const int m = c.sectionPolygon(p0, u, spx, spy, spz);
        const int64_t base = vOff(i);
        for (int k = 0; k < m; ++k) {
          verts(3 * (base + k)) = s[0] + spx[k];
          verts(3 * (base + k) + 1) = s[1] + spy[k];
          verts(3 * (base + k) + 2) = s[2] + spz[k];
        }
        const int p = polyOff(i);
        offset(p) = base;
        volume(p) = c.volumePerVertex();
        seed(p) = i;
      });
  Kokkos::fence();
  detail::statusCounts(status, N, res.numOverflow, res.numIncomplete);
  return res;
}
 
}  // namespace peclet::voro
 
#endif  // PECLET_VORO_PORE_CELLS_HPP