wall.hpp

Go to the documentation of this file.

 
#ifndef PECLET_VORO_ENERGY_WALL_HPP
#define PECLET_VORO_ENERGY_WALL_HPP
 
#include <Kokkos_Core.hpp>
#include <type_traits>
 
#include "peclet/voro/energy/route.hpp"
#include "peclet/voro/sdf.hpp"
 
namespace peclet::voro::energy {
 
template <class Real, class Policy = Voronoi, class Sdf = NoSdf>
Real wallEnergyForce(const TessellationView<Real>& view,
                     const Kokkos::View<int*, peclet::core::MemSpace>& types,
                     const Kokkos::View<Real*, peclet::core::MemSpace>& sigmaS, const Sdf& sdf,
                     const Kokkos::View<Real*, peclet::core::MemSpace>& pos,
                     const Kokkos::View<Real*, peclet::core::MemSpace>& weight, Real L,
                     const Kokkos::View<Real*, peclet::core::MemSpace>& force,
                     const Kokkos::View<Real*, peclet::core::MemSpace>& forceW = {}) {
  using Exec = peclet::core::ExecSpace;
  const int N = view.numCells();
  Real E = 0;
  Kokkos::parallel_reduce(
      "energy.wall", Kokkos::RangePolicy<Exec>(0, N),
      KOKKOS_LAMBDA(const int i, Real& acc) {
        const Real sig = sigmaS(types(i));
        if (sig == Real(0))
          return;
        Real gl[detail::kMaxLocalFacets * 3];
        const int b = view.facetBegin(i), nf = view.facetEnd(i) - b;
        if (nf > detail::kMaxLocalFacets)
          return;
        for (int k = 0; k < 3 * nf; ++k)
          gl[k] = Real(0);
        bool any = false;
        const bool fd = view.hasWallFD();  // exact wall-plane part (sdfWallFD) instead of the chain
        for (int f = b; f < b + nf; ++f) {
          if (view.facetNbr(f) != kBoundaryFacet)
            continue;
          any = true;
          const Real ax = view.area(f, 0), ay = view.area(f, 1), az = view.area(f, 2);
          acc += sig * Kokkos::sqrt(ax * ax + ay * ay + az * az);
          for (int e = view.edgeBegin(f); e < view.edgeEnd(f); ++e) {
            const int p = view.edgePartner(e) - b;
            if (fd && view.facetNbr(b + p) == kBoundaryFacet)
              continue;  // wall-plane dependence comes from the FD part below
            for (int c = 0; c < 3; ++c)
              gl[3 * p + c] += sig * view.areaGrad(e, c);
          }
        }
        if (any) {
          detail::routeCellGradient<Real, Policy>(view, i, gl, nf, pos, weight, L, sdf, force,
                                                  forceW);
          if (fd)
            for (int c = 0; c < 3; ++c)
              Kokkos::atomic_add(&force(3 * i + c), sig * view.wallDA(i, c));
        }
      },
      E);
  return E;
}
 
}  // namespace peclet::voro::energy
 
#endif  // PECLET_VORO_ENERGY_WALL_HPP