subset_gather.hpp

Go to the documentation of this file.

 
#ifndef PECLET_VORO_SUBSET_GATHER_HPP
#define PECLET_VORO_SUBSET_GATHER_HPP
 
#include <Kokkos_Core.hpp>
#include <string>
 
#include "peclet/core/common/view.hpp"
#include "peclet/voro/tess_grid.hpp"
#include "peclet/voro/tessellator.hpp"  // CellBuilder, StatusBit, NoSdf
 
namespace peclet::voro {
 
template <class Real>
struct SubsetGatherResult {
  Kokkos::View<int*, peclet::core::MemSpace>
      status;  // N : per-cell StatusBit mask (subset entries written)
};
 
template <class Real, bool Weighted = false, bool TrackAdj = false, class Sdf = NoSdf>
SubsetGatherResult<Real> subsetGather(
    const TessGrid<Real>& grid, const Kokkos::View<int*, peclet::core::MemSpace>& indices,
    int nSubset, const Kokkos::View<int*, peclet::core::MemSpace>& outNp,
    const Kokkos::View<int*, peclet::core::MemSpace>& outNt,
    const Kokkos::View<int*, peclet::core::MemSpace>& outPnbr,
    const Kokkos::View<unsigned*, peclet::core::MemSpace>& outTri,
    const Kokkos::View<Real*, peclet::core::MemSpace>& cellVol,
    const Kokkos::View<unsigned char*, peclet::core::MemSpace>& outPoke4 = {}, Sdf sdf = {},
    bool withForceGeom = false, WallStore<Real> outWall = {},
    Kokkos::View<int*, peclet::core::MemSpace> outNear = {},
    Kokkos::View<int*, peclet::core::MemSpace> outNearCnt = {}, int nearCap = 0,
    Real nearMargin = Real(0)) {
  using peclet::core::MemSpace;
  using Exec = peclet::core::ExecSpace;
  using Builder = CellBuilder<Real, Weighted, Sdf, TrackAdj>;
  const int N = grid.N;
  using Kokkos::view_alloc;
  using Kokkos::WithoutInitializing;
 
  // Power reach needs the global max weight (see CellBuilder::blockReachSq); 0 for Voronoi. The
  // grid carries the per-slot weights (wSorted), so reduce over those.
  Real wMaxAll = Real(0);
  if constexpr (Weighted) {
    if (grid.wSorted.extent(0) == static_cast<size_t>(N)) {
      auto ws = grid.wSorted;
      Kokkos::parallel_reduce(
          "subset.wmax", Kokkos::RangePolicy<Exec>(0, N),
          KOKKOS_LAMBDA(int idx, Real& m) { m = ws(idx) > m ? ws(idx) : m; },
          Kokkos::Max<Real>(wMaxAll));
    }
  }
 
  SubsetGatherResult<Real> res;
  res.status = Kokkos::View<int*, MemSpace>("subset.status", N);  // zero-init (unwritten = kOk)
 
  // Throwaway facet over-buffer: buildCell reserves a contiguous CSR range per cell and writes the
  // neighbour id (+ optional geometry) there. Sized for the subset (mean ~15.5 faces/cell +
  // headroom); the topology store is the published Phase-1 output, this buffer is just buildCell's
  // scratch.
  Kokkos::View<int*, MemSpace> facetCount("subset.facetCount", N);
  Kokkos::View<int*, MemSpace> cellFacetBase("subset.cellFacetBase", N);
  // Worst case per cell (MAXF_TMP): a subset of wall-clipped cells carries far more faces than the
  // Poisson mean, and the buffer is only scratch.
  const size_t facetCap = (size_t)(nSubset > 0 ? nSubset : 1) * (size_t)Builder::MAXF_TMP;
  Kokkos::View<int*, MemSpace> oNbr(view_alloc(std::string("subset.oNbr"), WithoutInitializing),
                                    facetCap);
  Kokkos::View<Real*, MemSpace> oArea(view_alloc(std::string("subset.oArea"), WithoutInitializing),
                                      facetCap * 3);
  Kokkos::View<Real*, MemSpace> oDV(view_alloc(std::string("subset.oDV"), WithoutInitializing),
                                    facetCap * 3);
  Kokkos::View<Real*, MemSpace> oConn(view_alloc(std::string("subset.oConn"), WithoutInitializing),
                                      facetCap * 3);
  Kokkos::View<int*, MemSpace> facetCursor("subset.facetCursor", 1);
 
  // Empty candidate-list outputs (skin emission is not wanted for the subset gather).
  Kokkos::View<int*, MemSpace> noCand, noCandCnt;
 
  Builder op{grid.binned,
             grid.posSorted,
             grid.wSorted,  // per-slot weights (empty for Voronoi; read only when Weighted)
             grid.gidSorted,
             grid.cellStart,
             grid.wlOff,
             grid.wlRmin,
             res.status,
             cellVol,
             facetCount,
             cellFacetBase,
             oNbr,
             oArea,
             oDV,
             oConn,
             facetCursor,
             grid.icx,
             grid.icy,
             grid.icz,
             grid.Lx,
             grid.Ly,
             grid.Lz,
             grid.minCsz,
             wMaxAll,
             grid.dimx,
             grid.dimy,
             grid.dimz,
             grid.sw,
             grid.nOff,
             grid.wlS,
             grid.useMorton,
             grid.haveGid,
             withForceGeom,
             facetCap,
             sdf,
             outNp,
             outNt,
             outPnbr,
             outTri,
             outPoke4,
             noCand,
             noCandCnt,
             /*emitTopo=*/true,
             /*emitCand=*/false,
             /*candCap=*/0,
             outWall,
             /*oEdgeOff=*/{},
             /*oEdgeCnt=*/{},
             /*oEdgeFacet=*/{},
             /*oEdgeGrad=*/{},
             /*oEdgeLen=*/{},
             /*edgeCursor=*/{},
             /*edgeCap=*/0,
             /*withAreaGrad=*/false,
             outNear,
             outNearCnt,
             nearCap,
             nearMargin,
             /*oWallDV=*/{},
             /*oWallDA=*/{},
             /*oMoment2=*/{}};
 
  auto slotOf = grid.slotOf;
  auto idx = indices;
  Kokkos::parallel_for(
      "subset.build", Kokkos::RangePolicy<Exec>(0, nSubset), KOKKOS_LAMBDA(const int s) {
        const int i = idx(s);  // original seed index to rebuild
        op.buildCell(
            slotOf(i));  // grid-sorted slot of that seed; writes outputs at binned(slot)==i
      });
  Kokkos::fence();
  return res;
}
 
}  // namespace peclet::voro
 
#endif  // PECLET_VORO_SUBSET_GATHER_HPP