transpose.hpp

Go to the documentation of this file.

 
#ifndef PECLET_VORO_TRANSPOSE_HPP
#define PECLET_VORO_TRANSPOSE_HPP
 
#include <Kokkos_Core.hpp>
 
#include "peclet/core/common/view.hpp"
#include "peclet/voro/tessellation_view.hpp"
 
namespace peclet::voro {
 
template <class Real>
struct AuxMaps {
  Kokkos::View<int*, peclet::core::MemSpace> recip;        // nFacets: reciprocal facet, or -1
  Kokkos::View<int*, peclet::core::MemSpace> cellOfFacet;  // nFacets: owning cell index
};
 
template <class Real>
AuxMaps<Real> buildAuxMaps(const TessellationView<Real>& view) {
  using Exec = peclet::core::ExecSpace;
  const int N = view.numCells();
  const int nF = view.numFacets();
  AuxMaps<Real> aux;
  aux.recip = Kokkos::View<int*, peclet::core::MemSpace>("recip", nF);
  aux.cellOfFacet = Kokkos::View<int*, peclet::core::MemSpace>("cellOfFacet", nF);
  auto recip = aux.recip;
  auto cellOfFacet = aux.cellOfFacet;
 
  Kokkos::parallel_for(
      "aux.cellOfFacet", Kokkos::RangePolicy<Exec>(0, N), KOKKOS_LAMBDA(const int i) {
        for (int f = view.facetBegin(i); f < view.facetEnd(i); ++f)
          cellOfFacet(f) = i;
      });
 
  Kokkos::parallel_for(
      "aux.recip", Kokkos::RangePolicy<Exec>(0, N), KOKKOS_LAMBDA(const int i) {
        const gid_t idi = view.cellSeed(i);
        for (int g = view.facetBegin(i); g < view.facetEnd(i); ++g) {
          const gid_t j = view.facetNbr(g);
          recip(g) = -1;
          if (j < 0 || j >= N)
            continue;
          int best = -1;
          Real bestErr = Real(1e300);
          for (int h = view.facetBegin((int)j); h < view.facetEnd((int)j); ++h) {
            if (view.facetNbr(h) != idi)
              continue;
            Real e = 0;  // area-vector negation picks the matching interface
            for (int c = 0; c < 3; ++c) {
              Real s = view.area(g, c) + view.area(h, c);
              e += s * s;
            }
            if (e < bestErr) {
              bestErr = e;
              best = h;
            }
          }
          recip(g) = best;
        }
      });
  Kokkos::fence();
  return aux;
}
 
}  // namespace peclet::voro
 
#endif  // PECLET_VORO_TRANSPOSE_HPP