voronoi_halo.hpp

Go to the documentation of this file.

 
#ifndef PECLET_VORO_MPI_VORONOI_HALO_HPP
#define PECLET_VORO_MPI_VORONOI_HALO_HPP
 
#include <array>
#include <cstddef>
#include <vector>
 
#include "peclet/core/common/mpi.hpp"
#include "peclet/core/decomp/block_decomposer.hpp"
#include "peclet/core/halo/particle_halo_topology.hpp"
#include "peclet/core/halo/particle_migrator.hpp"
 
namespace peclet::voro {
namespace mpi {
 
template <class Real>
class VoronoiHalo {
 public:
  using Vec3 = std::array<Real, 3>;
 
  struct Gathered {
    std::vector<Vec3> pos;     
    std::vector<long> gid;     
    std::vector<Real> weight;  
    int nOwned = 0;
  };
 
  void init(std::array<Real, 3> origin, std::array<Real, 3> size, std::array<long, 3> gsize,
            std::array<bool, 3> periodic, MPI_Comm comm = MPI_COMM_WORLD) {
    comm_ = comm;
    int sz = 1;
    MPI_Comm_rank(comm_, &rank_);
    MPI_Comm_size(comm_, &sz);
    dec_.init(static_cast<std::size_t>(sz), peclet::core::IVec<3>{gsize[0], gsize[1], gsize[2]});
    peclet::core::halo::DomainMap<3> map;
    for (int i = 0; i < 3; ++i) {
      map.origin[i] = origin[i];
      map.cellSize[i] = size[i] / static_cast<double>(gsize[i]);
      map.periodic[i] = periodic[i];
    }
    mig_.init(dec_, rank_, map, comm_);
    halo_.init(mig_);
  }
 
  int rank() const { return rank_; }
  int size() const {
    int s = 1;
    MPI_Comm_size(comm_, &s);
    return s;
  }
 
  int ownerOf(const Vec3& x) const {
    peclet::core::Vec<3> v{x[0], x[1], x[2]};
    return mig_.ownerOf(v);
  }
 
  Gathered gather(const std::vector<Vec3>& ownedPos, const std::vector<long>& ownedGid,
                  const std::vector<Real>& ownedWeight, double rcut) {
    const int nOwned = static_cast<int>(ownedPos.size());
    std::vector<peclet::core::Vec<3>> pv(nOwned);
    for (int i = 0; i < nOwned; ++i)
      pv[i] = peclet::core::Vec<3>{ownedPos[i][0], ownedPos[i][1], ownedPos[i][2]};
    // includePeriodicSelf=false: the tessellator is periodic-native (minimal image),
    // so two same-rank seeds interact across a periodic face directly — only
    // cross-rank neighbours (which the build still gathers, periodic images
    // included) need to travel as ghosts. Self-ghosts would just duplicate owned
    // seeds (wrapping back onto them) and inflate the set.
    halo_.build(pv, rcut, /*includePeriodicSelf=*/false);
    const int ng = static_cast<int>(halo_.numGhost());
 
    // Forward owner global-id and weight onto the ghost copies.
    std::vector<long> ghostGid(ng);
    std::vector<Real> ghostW(ng);
    halo_.forward(ownedGid.data(), ghostGid.data());
    halo_.forward(ownedWeight.data(), ghostW.data());
    const std::vector<peclet::core::Vec<3>>& gpos = halo_.ghostPositions();
 
    Gathered g;
    g.nOwned = nOwned;
    g.pos.resize(nOwned + ng);
    g.gid.resize(nOwned + ng);
    g.weight.resize(nOwned + ng);
    for (int i = 0; i < nOwned; ++i) {
      g.pos[i] = ownedPos[i];
      g.gid[i] = ownedGid[i];
      g.weight[i] = ownedWeight[i];
    }
    for (int j = 0; j < ng; ++j) {
      g.pos[nOwned + j] = Vec3{gpos[j][0], gpos[j][1], gpos[j][2]};
      g.gid[nOwned + j] = ghostGid[j];
      g.weight[nOwned + j] = ghostW[j];
    }
    nGhost_ = ng;
    return g;
  }
 
  int numGhost() const { return nGhost_; }
 
  void refreshPositions(const std::vector<Vec3>& ownedPos, std::vector<Vec3>& out) {
    const int nOwned = static_cast<int>(ownedPos.size());
    const int ng = nGhost_;
    std::vector<Real> ox(nOwned), oy(nOwned), oz(nOwned), gx(ng), gy(ng), gz(ng);
    for (int i = 0; i < nOwned; ++i) {
      ox[i] = ownedPos[i][0];
      oy[i] = ownedPos[i][1];
      oz[i] = ownedPos[i][2];
    }
    halo_.forward(ox.data(), gx.data());
    halo_.forward(oy.data(), gy.data());
    halo_.forward(oz.data(), gz.data());
    out.resize(nOwned + ng);
    for (int i = 0; i < nOwned; ++i)
      out[i] = ownedPos[i];
    for (int j = 0; j < ng; ++j)
      out[nOwned + j] = Vec3{gx[j], gy[j], gz[j]};
  }
 
 private:
  int nGhost_ = 0;
  MPI_Comm comm_ = MPI_COMM_WORLD;
  int rank_ = 0;
  peclet::core::decomp::BlockDecomposer<3> dec_;
  peclet::core::halo::ParticleMigrator<3> mig_;
  peclet::core::halo::ParticleHaloTopology<3> halo_;
};
 
}  // namespace mpi
}  // namespace peclet::voro
 
#endif  // PECLET_VORO_MPI_VORONOI_HALO_HPP