interface.hpp

Go to the documentation of this file.

 
#ifndef PECLET_VORO_PHYSICS_INTERFACE_HPP
#define PECLET_VORO_PHYSICS_INTERFACE_HPP
 
#include <Kokkos_Core.hpp>
 
#include "peclet/core/common/view.hpp"
#include "peclet/voro/tessellation_view.hpp"
 
namespace peclet::voro {
namespace physics {
 
template <class Real>
Real interfaceEnergy(const TessellationView<Real>& view,
                     const Kokkos::View<int*, peclet::core::MemSpace>& types,
                     const Kokkos::View<Real*, peclet::core::MemSpace>& tension, int nTypes) {
  using Exec = peclet::core::ExecSpace;
  const int N = view.numCells();
  Real e = 0;
  Kokkos::parallel_reduce(
      "intf.energy", Kokkos::RangePolicy<Exec>(0, N),
      KOKKOS_LAMBDA(const int i, Real& acc) {
        const int ti = types(i);
        if (ti == 0)
          return;
        for (int f = view.facetBegin(i); f < view.facetEnd(i); ++f) {
          const int j = view.facetNbr(f);
          if (j < 0 || j >= N)
            continue;
          const int tj = types(j);
          if (tj < ti) {
            Real ax = view.area(f, 0), ay = view.area(f, 1), az = view.area(f, 2);
            acc += tension(ti * nTypes + tj) * Kokkos::sqrt(ax * ax + ay * ay + az * az);
          }
        }
      },
      e);
  return e;
}
 
}  // namespace physics
}  // namespace peclet::voro
 
#endif  // PECLET_VORO_PHYSICS_INTERFACE_HPP