interface.cpp

Go to the documentation of this file.

//
// This file is subject to the terms and conditions defined in
// file 'LICENSE', which is part of this source code package.
//
 
// Define before including interpolate.hpp so the extern-template suppression
// block is skipped and this TU can provide explicit instantiations.
#define PACMAN_FE_ETI_COMPILATION
 
#include "interface.hpp"
 
#include <Kokkos_Core.hpp>
#include <Kokkos_Profiling_ScopedRegion.hpp>
#include <stdexcept>
#include <string>
#include <unordered_map>
 
#include "common/concepts.hpp"
#include "common/transfer.hxx"
#include "interpolate.hpp"
 
// ---------------------------------------------------------------------------
// Internal helpers (anonymous namespace)
// ---------------------------------------------------------------------------
 
namespace {
 
template <typename ExecSpace, PACMAN::int_t Dim>
void
RunInterpolate(PACMAN::method_t method, PACMAN::coordinates_t *pSourcePoints,
               PACMAN::int_t nSourcePoints, PACMAN::fp_t *pSourceValues,
               PACMAN::int_t *pConnVal, PACMAN::int_t connValSize,
               PACMAN::offset_t *pConnOff, PACMAN::int_t connOffSize,
               PACMAN::cell_t *pCellTypes, PACMAN::coordinates_t *pTargetPoints,
               PACMAN::int_t nTargetPoints, bool fortranIndexing,
               PACMAN::int_t nComponents, PACMAN::fp_t *pTargetValues,
               PACMAN::int_t *pTargetStatus) {
  PACMAN::Transfer<ExecSpace, Dim> transfer(method);
  PACMAN::SetupTransferClass(transfer, nSourcePoints, connValSize, connOffSize,
                             nTargetPoints, nComponents, pSourcePoints,
                             pSourceValues,
                             pConnVal, pConnOff, pCellTypes, pTargetPoints, fortranIndexing);
  PACMAN::Interpolate(transfer);
 
  auto unmanaged_host_ts = Kokkos::View<
      PACMAN::TransferStatus *, Kokkos::DefaultHostExecutionSpace::memory_space,
      Kokkos::MemoryTraits<Kokkos::Unmanaged>>(
      reinterpret_cast<PACMAN::TransferStatus *>(pTargetStatus), nTargetPoints);
 
  PACMAN::CopyToRowMajorHost(transfer.targetValues, pTargetValues);
  Kokkos::deep_copy(unmanaged_host_ts, transfer.targetStatus);
}
 
template <PACMAN::int_t Dim>
requires PACMAN::IsValidDim<Dim> void
Dispatch(unsigned char execSpace, PACMAN::method_t method,
         PACMAN::coordinates_t *pSourcePoints, PACMAN::int_t nSourcePoints,
         PACMAN::fp_t *pSourceValues, PACMAN::int_t *pConnVal,
         PACMAN::int_t connValSize, PACMAN::offset_t *pConnOff,
         PACMAN::int_t connOffSize, PACMAN::cell_t *pCellTypes,
         PACMAN::coordinates_t *pTargetPoints, PACMAN::int_t nTargetPoints,
         bool fortranIndexing, PACMAN::int_t nComponents,
         PACMAN::fp_t *pTargetValues, PACMAN::int_t *pTargetStatus) {
      std::visit(
      [&](auto execSpaceObj) {
        RunInterpolate<decltype(execSpaceObj), Dim>(
            method, pSourcePoints, nSourcePoints, pSourceValues, pConnVal,
            connValSize, pConnOff, connOffSize, pCellTypes, pTargetPoints,
          nTargetPoints, fortranIndexing, nComponents, pTargetValues,
          pTargetStatus);
      },
      PACMAN::MakeExecSpace(execSpace));
}
 
} // anonymous namespace
 
// ---------------------------------------------------------------------------
// Public API implementation
// ---------------------------------------------------------------------------
 
namespace PACMAN {
 
void
fe_interpolate(int_t spaceDimension, unsigned char execSpace, method_t method,
               coordinates_t *sourcePoints, int_t nSourcePoints,
               fp_t *sourceValues, int_t *connVal, int_t connValSize,
               offset_t *connOff, int_t connOffSize, cell_t *cellTypes,
               coordinates_t *targetPoints, int_t nTargetPoints,
               fp_t *targetValues, int_t *targetStatus, bool fortranIndexing,
               int_t nComponents) {
  const std::string _region_name = "PACMAN::fe_interpolate";
  const Kokkos::Profiling::ScopedRegion region(_region_name);
 
  // LCOV_EXCL_START
  if (sourcePoints == nullptr && nSourcePoints > 0)
    throw std::invalid_argument("sourcePoints is null but nSourcePoints > 0");
  if (targetPoints == nullptr && nTargetPoints > 0)
    throw std::invalid_argument("targetPoints is null but nTargetPoints > 0");
  if (sourceValues == nullptr && nSourcePoints > 0)
    throw std::invalid_argument("sourceValues is null but nSourcePoints > 0");
  if (targetValues == nullptr && nTargetPoints > 0)
    throw std::invalid_argument("targetValues is null but nTargetPoints > 0");
  if (targetStatus == nullptr && nTargetPoints > 0)
    throw std::invalid_argument("targetStatus is null but nTargetPoints > 0");
  if (nComponents < 1)
    throw std::invalid_argument("nComponents must be greater than zero");
  // LCOV_EXCL_STOP
 
  switch (spaceDimension) {
  case 1:
    return Dispatch<1>(execSpace, method, sourcePoints, nSourcePoints,
                       sourceValues, connVal, connValSize, connOff, connOffSize,
                       cellTypes, targetPoints, nTargetPoints, fortranIndexing,
                       nComponents, targetValues, targetStatus);
  case 2:
    return Dispatch<2>(execSpace, method, sourcePoints, nSourcePoints,
                       sourceValues, connVal, connValSize, connOff, connOffSize,
                       cellTypes, targetPoints, nTargetPoints, fortranIndexing,
                       nComponents, targetValues, targetStatus);
  case 3:
    return Dispatch<3>(execSpace, method, sourcePoints, nSourcePoints,
                       sourceValues, connVal, connValSize, connOff, connOffSize,
                       cellTypes, targetPoints, nTargetPoints, fortranIndexing,
                       nComponents, targetValues, targetStatus);
  // LCOV_EXCL_START
  default:
    throw std::runtime_error(
        "The dimension of the points can only be: 1, 2 or 3.\n");
  // LCOV_EXCL_STOP
  }
}
 
void vtk_to_pacman_cell_type(const int_t *vtkTypes, cell_t *pacmanTypes,
                             int_t n) {
  static const std::unordered_map<int_t, CellType> vtk_to_pacman = {
      // 0D
      {1, CellType::VTK_VERTEX},
      // 1D
      {3, CellType::VTK_LINE},
      {21, CellType::VTK_QUADRATIC_EDGE},
      // 2D
      {5, CellType::VTK_TRIANGLE},
      {9, CellType::VTK_QUAD},
      {22, CellType::VTK_QUADRATIC_TRIANGLE},
      {23, CellType::VTK_QUADRATIC_QUAD},
      // 3D
      {10, CellType::VTK_TETRA},
      {12, CellType::VTK_HEXAHEDRON},
      {13, CellType::VTK_WEDGE},
      {14, CellType::VTK_PYRAMID},
      {24, CellType::VTK_QUADRATIC_TETRA},
      {25, CellType::VTK_QUADRATIC_HEXAHEDRON},
      {26, CellType::VTK_QUADRATIC_WEDGE},
      {27, CellType::VTK_QUADRATIC_PYRAMID},
  };
  for (int_t i = 0; i < n; ++i) {
    auto it = vtk_to_pacman.find(vtkTypes[i]);
    // LCOV_EXCL_START
    if (it == vtk_to_pacman.end()) {
      throw std::runtime_error("Unsupported VTK cell type: " +
                               std::to_string(vtkTypes[i]));
    }
    // LCOV_EXCL_STOP
    pacmanTypes[i] = static_cast<cell_t>(it->second);
  }
}
 
} // namespace PACMAN