tessellation_view.hpp

Go to the documentation of this file.

 
#ifndef PECLET_VORO_TESSELLATION_VIEW_HPP
#define PECLET_VORO_TESSELLATION_VIEW_HPP
 
#include <cmath>
#include <cstdint>
#include <Kokkos_Core.hpp>
#include <map>
#include <vector>
 
#include "peclet/core/common/view.hpp"
namespace peclet::voro {
 
using gid_t = std::int32_t;
 
inline gid_t toGid(std::uint32_t id) {
  return static_cast<gid_t>(static_cast<std::int32_t>(id));
}
 
template <class Real>
struct HostTessellation {
  int nCells = 0;
  int nFacets = 0;
  std::vector<int> cellFacetOffset;  
  std::vector<gid_t> cellSeedId;     
  std::vector<Real> cellVolume;      
  std::vector<gid_t> facetNeighbor;  
  std::vector<Real> facetArea;       
  std::vector<Real> facetConnect;    
  std::vector<Real> facetConnVec;    
};
 
template <class Real>
std::vector<int> buildReciprocalMap(const HostTessellation<Real>& t) {
  std::map<gid_t, int> seedToCell;
  for (int i = 0; i < t.nCells; ++i)
    seedToCell[t.cellSeedId[i]] = i;
  std::vector<int> recip(t.nFacets, -1);
  for (int i = 0; i < t.nCells; ++i) {
    const gid_t idi = t.cellSeedId[i];
    for (int g = t.cellFacetOffset[i]; g < t.cellFacetOffset[i + 1]; ++g) {
      const gid_t j = t.facetNeighbor[g];
      if (j < 0)
        continue;  // boundary
      auto it = seedToCell.find(j);
      if (it == seedToCell.end())
        continue;
      const int cj = it->second;
      int best = -1;
      Real bestErr = Real(1e300);
      for (int h = t.cellFacetOffset[cj]; h < t.cellFacetOffset[cj + 1]; ++h) {
        if (t.facetNeighbor[h] != idi)
          continue;
        Real e = 0;  // prefer the facet whose area best negates g's (same interface)
        for (int c = 0; c < 3; ++c) {
          Real s = t.facetArea[3 * g + c] + t.facetArea[3 * h + c];
          e += s * s;
        }
        if (e < bestErr) {
          bestErr = e;
          best = h;
        }
      }
      recip[g] = best;
    }
  }
  return recip;
}
 
template <class Real>
struct TessellationView {
  // facetBegin(i) = cellFacetOffset(i); facetEnd(i) = cellFacetOffset(i)+cellFacetCount(i).
  // The device tessellator packs the CSR in finish-order via an atomic cursor (so
  // cellFacetOffset is a per-cell base, not a prefix sum), hence the explicit count;
  // the legacy upload path fills count from its prefix-sum offsets. Consumers use
  // facetBegin/facetEnd and are agnostic to which.
  Kokkos::View<int*, peclet::core::MemSpace> cellFacetOffset;  // nCells (per-cell facet base)
  Kokkos::View<int*, peclet::core::MemSpace> cellFacetCount;   // nCells
  Kokkos::View<gid_t*, peclet::core::MemSpace> cellSeedId;     // nCells
  Kokkos::View<Real*, peclet::core::MemSpace> cellVolume;      // nCells
  Kokkos::View<gid_t*, peclet::core::MemSpace> facetNeighbor;  // nFacets
  Kokkos::View<Real*, peclet::core::MemSpace> facetArea;       // 3*nFacets
  Kokkos::View<Real*, peclet::core::MemSpace> facetConnect;    // 3*nFacets (dV)
  Kokkos::View<Real*, peclet::core::MemSpace> facetConnVec;    // 3*nFacets
  // Vertex CSR (seed-relative positions) — published when a consumer needs cell
  // geometry (rendering, the interface-tension force, future Hamiltonian solver);
  // empty otherwise.
  Kokkos::View<int*, peclet::core::MemSpace> cellVertexOffset;  // nCells+1 (empty if unused)
  Kokkos::View<Real*, peclet::core::MemSpace> vertexPos;        // 3*nVerts
 
  KOKKOS_INLINE_FUNCTION int numCells() const { return static_cast<int>(cellSeedId.extent(0)); }
  KOKKOS_INLINE_FUNCTION bool hasVertices() const { return cellVertexOffset.extent(0) > 0; }
  KOKKOS_INLINE_FUNCTION int vertexBegin(int i) const { return cellVertexOffset(i); }
  KOKKOS_INLINE_FUNCTION int vertexEnd(int i) const { return cellVertexOffset(i + 1); }
  KOKKOS_INLINE_FUNCTION Real vertex(int v, int c) const { return vertexPos(3 * v + c); }
  KOKKOS_INLINE_FUNCTION int numFacets() const { return static_cast<int>(facetNeighbor.extent(0)); }
  KOKKOS_INLINE_FUNCTION gid_t cellSeed(int i) const { return cellSeedId(i); }
  KOKKOS_INLINE_FUNCTION Real volume(int i) const { return cellVolume(i); }
  KOKKOS_INLINE_FUNCTION int facetBegin(int i) const { return cellFacetOffset(i); }
  KOKKOS_INLINE_FUNCTION int facetEnd(int i) const {
    return cellFacetOffset(i) + cellFacetCount(i);
  }
  KOKKOS_INLINE_FUNCTION gid_t facetNbr(int f) const { return facetNeighbor(f); }
  KOKKOS_INLINE_FUNCTION Real area(int f, int c) const { return facetArea(3 * f + c); }
  KOKKOS_INLINE_FUNCTION Real connect(int f, int c) const { return facetConnect(3 * f + c); }
  KOKKOS_INLINE_FUNCTION Real connVec(int f, int c) const { return facetConnVec(3 * f + c); }
};
 
namespace detail {
template <class T>
Kokkos::View<T*, peclet::core::MemSpace> viewFrom(const std::vector<T>& h,
                                                  const std::string& label) {
  // NOTE: pass a std::string label — a bare const char* is interpreted by
  // Kokkos::view_alloc as a user pointer-to-memory, not an allocation label.
  Kokkos::View<T*, peclet::core::MemSpace> d(Kokkos::view_alloc(label, Kokkos::WithoutInitializing),
                                             h.size());
  if (!h.empty()) {
    using Host =
        Kokkos::View<const T*, Kokkos::HostSpace, Kokkos::MemoryTraits<Kokkos::Unmanaged> >;
    Host hv(h.data(), h.size());
    Kokkos::deep_copy(d, hv);
  }
  return d;
}
}  // namespace detail
 
template <class Real>
TessellationView<Real> upload(const HostTessellation<Real>& h) {
  TessellationView<Real> v;
  // Legacy CSR is a prefix sum (nCells+1); derive the per-cell count so facetEnd
  // works (the device path instead writes base + count directly).
  v.cellFacetOffset = detail::viewFrom(h.cellFacetOffset, "tess.cellFacetOffset");
  {
    const int nC = static_cast<int>(h.cellSeedId.size());
    std::vector<int> cnt(nC);
    for (int i = 0; i < nC; ++i)
      cnt[i] = h.cellFacetOffset[i + 1] - h.cellFacetOffset[i];
    v.cellFacetCount = detail::viewFrom(cnt, "tess.cellFacetCount");
  }
  v.cellSeedId = detail::viewFrom(h.cellSeedId, "tess.cellSeedId");
  v.cellVolume = detail::viewFrom(h.cellVolume, "tess.cellVolume");
  v.facetNeighbor = detail::viewFrom(h.facetNeighbor, "tess.facetNeighbor");
  v.facetArea = detail::viewFrom(h.facetArea, "tess.facetArea");
  v.facetConnect = detail::viewFrom(h.facetConnect, "tess.facetConnect");
  v.facetConnVec = detail::viewFrom(h.facetConnVec, "tess.facetConnVec");
  return v;
}
 
}  // namespace peclet::voro
 
#endif  // PECLET_VORO_TESSELLATION_VIEW_HPP