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
  // Facet-edge AREA-JACOBIAN CSR (Voronoi methods plan, rung A3; opt-in `withAreaGrad`): for
  // facet f, entries e in [edgeBegin(f), edgeEnd(f)) hold ∂A_f/∂n_l for every plane l the
  // facet's area depends on — slot 0 is the facet's OWN plane (l = f, the diagonal block), the
  // rest are the facets sharing an edge with f (each edge is where two of the cell's planes
  // meet; a face's area depends on its own plane and on its edge-neighbours' planes only).
  // `edgeFacet(e)` is the global facet index of plane l's facet in the SAME cell (its neighbour
  // id / wall sentinel is facetNbr(edgeFacet(e))); `areaGrad(e, c)` the 3-vector ∂A_f/∂n_l.
  // A_f is the facet's area MAGNITUDE (|facetArea|). Empty unless requested.
  Kokkos::View<int*, peclet::core::MemSpace> facetEdgeOffset;  // nFacets
  Kokkos::View<int*, peclet::core::MemSpace> facetEdgeCount;   // nFacets
  Kokkos::View<int*, peclet::core::MemSpace> edgeFacet;        // nEdges
  Kokkos::View<Real*, peclet::core::MemSpace> edgeAreaGrad;    // 3*nEdges
  // Track C2a′: the Voronoi EDGE length of each facet-edge slot (0 for the self slot): the edge
  // where the facet's plane meets the partner facet's plane — the dual of the Delaunay triangle
  // (seed, facet neighbour, partner neighbour); the DEC curl-curl viscous term's Hodge star.
  Kokkos::View<Real*, peclet::core::MemSpace> edgeLength;  // nEdges
 
  // Rung A1 (force half; opt-in `withWallFD`): per cell, the WALL part of dV/dx and dA_wall/dx
  // by in-kernel central differences of the SDF clip (sdfWallFD) — exact for the clip actually
  // used (sagitta placement included), where the seed-foot chain is the tangent-plane model.
  // Zero for cells without wall planes. Empty unless requested.
  Kokkos::View<Real*, peclet::core::MemSpace> cellWallDV;  // 3*nCells
  Kokkos::View<Real*, peclet::core::MemSpace> cellWallDA;  // 3*nCells
  KOKKOS_INLINE_FUNCTION bool hasWallFD() const { return cellWallDV.extent(0) > 0; }
  KOKKOS_INLINE_FUNCTION Real wallDV(int i, int c) const { return cellWallDV(3 * i + c); }
  KOKKOS_INLINE_FUNCTION Real wallDA(int i, int c) const { return cellWallDA(3 * i + c); }
 
  // Track B, rung B1 (opt-in `withMoments`): per facet ∫_face |s|² dA, the second moment of the
  // face polygon about the seed's foot point on its plane — what the Lloyd (centroidal) energy
  // value needs (a pyramid's ∫|y|² dy = h (h² A + M2)/5). Empty unless requested.
  Kokkos::View<Real*, peclet::core::MemSpace> facetMoment2;  // nFacets
  KOKKOS_INLINE_FUNCTION bool hasMoments() const { return facetMoment2.extent(0) > 0; }
  KOKKOS_INLINE_FUNCTION Real moment2(int f) const { return facetMoment2(f); }
 
  KOKKOS_INLINE_FUNCTION bool hasAreaGrad() const { return facetEdgeOffset.extent(0) > 0; }
  KOKKOS_INLINE_FUNCTION int numEdges() const { return static_cast<int>(edgeFacet.extent(0)); }
  KOKKOS_INLINE_FUNCTION int edgeBegin(int f) const { return facetEdgeOffset(f); }
  KOKKOS_INLINE_FUNCTION int edgeEnd(int f) const { return facetEdgeOffset(f) + facetEdgeCount(f); }
  KOKKOS_INLINE_FUNCTION int edgePartner(int e) const { return edgeFacet(e); }
  KOKKOS_INLINE_FUNCTION Real areaGrad(int e, int c) const { return edgeAreaGrad(3 * e + c); }
  KOKKOS_INLINE_FUNCTION Real edgeLen(int e) const { return edgeLength(e); }
  KOKKOS_INLINE_FUNCTION bool hasEdgeLength() const { return edgeLength.extent(0) > 0; }
 
  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