lloyd.hpp

Go to the documentation of this file.

 
#ifndef PECLET_VORO_ENERGY_LLOYD_HPP
#define PECLET_VORO_ENERGY_LLOYD_HPP
 
#include <Kokkos_Core.hpp>
 
#include "peclet/voro/energy/route.hpp"
#include "peclet/voro/tessellation_view.hpp"
 
namespace peclet::voro::energy {
 
template <class Real>
void cellCentroids(const TessellationView<Real>& view,
                   const Kokkos::View<Real*, peclet::core::MemSpace>& pos,
                   const Kokkos::View<Real*, peclet::core::MemSpace>& out) {
  using Exec = peclet::core::ExecSpace;
  const int N = view.numCells();
  Kokkos::parallel_for(
      "energy.centroid", Kokkos::RangePolicy<Exec>(0, N), KOKKOS_LAMBDA(const int i) {
        Real m[3] = {0, 0, 0}, V = 0;
        for (int f = view.facetBegin(i); f < view.facetEnd(i); ++f) {
          const Real ax = view.area(f, 0), ay = view.area(f, 1), az = view.area(f, 2);
          const Real A = Kokkos::sqrt(ax * ax + ay * ay + az * az);
          if (!(A > Real(0)))
            continue;
          const Real rx = view.connVec(f, 0), ry = view.connVec(f, 1), rz = view.connVec(f, 2);
          const Real rl = Kokkos::sqrt(rx * rx + ry * ry + rz * rz);
          const Real h = Real(0.5) * rl;  // foot distance (conn = 2 n)
          const Real s = rl / A;
          const Real cx = rx - view.connect(f, 0) * s, cy = ry - view.connect(f, 1) * s,
                     cz = rz - view.connect(f, 2) * s;  // face centroid, seed-relative
          const Real Vf = h * A / Real(3);
          V += Vf;
          m[0] += Vf * Real(0.75) * cx;
          m[1] += Vf * Real(0.75) * cy;
          m[2] += Vf * Real(0.75) * cz;
        }
        const Real iv = V > Real(0) ? Real(1) / V : Real(0);
        for (int c = 0; c < 3; ++c)
          out(3 * i + c) = pos(3 * i + c) + m[c] * iv;
      });
}
 
template <class Real>
Real lloydEnergyForce(const TessellationView<Real>& view,
                      [[maybe_unused]] const Kokkos::View<Real*, peclet::core::MemSpace>& pos,
                      Real gamma, const Kokkos::View<Real*, peclet::core::MemSpace>& force) {
  using Exec = peclet::core::ExecSpace;
  const int N = view.numCells();
  const bool haveM2 = view.hasMoments();
  Real E = 0;
  Kokkos::parallel_reduce(
      "energy.lloyd", Kokkos::RangePolicy<Exec>(0, N),
      KOKKOS_LAMBDA(const int i, Real& acc) {
        Real m[3] = {0, 0, 0}, V = 0, e = 0;
        for (int f = view.facetBegin(i); f < view.facetEnd(i); ++f) {
          const Real ax = view.area(f, 0), ay = view.area(f, 1), az = view.area(f, 2);
          const Real A = Kokkos::sqrt(ax * ax + ay * ay + az * az);
          if (!(A > Real(0)))
            continue;
          const Real rx = view.connVec(f, 0), ry = view.connVec(f, 1), rz = view.connVec(f, 2);
          const Real rl = Kokkos::sqrt(rx * rx + ry * ry + rz * rz);
          const Real h = Real(0.5) * rl;
          const Real s = rl / A;
          const Real cx = rx - view.connect(f, 0) * s, cy = ry - view.connect(f, 1) * s,
                     cz = rz - view.connect(f, 2) * s;
          const Real Vf = h * A / Real(3);
          V += Vf;
          m[0] += Vf * Real(0.75) * cx;
          m[1] += Vf * Real(0.75) * cy;
          m[2] += Vf * Real(0.75) * cz;
          if (haveM2)
            e += h * (h * h * A + view.moment2(f)) / Real(5);
        }
        // gradient: 2 V (x − c) = −2 m  (m = V (c − x))
        for (int c = 0; c < 3; ++c)
          Kokkos::atomic_add(&force(3 * i + c), Real(-2) * gamma * m[c]);
        acc += gamma * e;
      },
      E);
  return haveM2 ? E : Real(-1);
}
 
}  // namespace peclet::voro::energy
 
#endif  // PECLET_VORO_ENERGY_LLOYD_HPP