kernel.h

// Copyright (C) 2025 Paul T. Kühner
//
// This file is part of DOLFINx (https://www.fenicsproject.org)
//
// SPDX-License-Identifier:    LGPL-3.0-or-later
 
#pragma once
 
#include <dolfinx/common/defines.h>
#include <dolfinx/common/types.h>
#include <ufcx.h>
 
namespace dolfinx::fem::impl
{
template <dolfinx::scalar T, std::floating_point U = scalar_value_t<T>>
using kernelptr_t = void (*)(T*, const T*, const T*, const U*, const int*,
                             const std::uint8_t*, void*);

template <dolfinx::scalar T, std::floating_point U = scalar_value_t<T>>
using kernel_t = std::function<void(T*, const T*, const T*, const U*,
                                    const int*, const std::uint8_t*, void*)>;

template <dolfinx::scalar T, std::floating_point U = scalar_value_t<T>>
constexpr kernel_t<T, U> extract_kernel(const ufcx_integral* integral)
{
  static_assert(std::is_same_v<U, scalar_value_t<T>>,
                "UFCx kernels require geometry type U == scalar_value_t<T>.");
 
  if constexpr (std::is_same_v<T, float>)
    return integral->tabulate_tensor_float32;
  else if constexpr (std::is_same_v<T, double>)
    return integral->tabulate_tensor_float64;
#ifndef DOLFINX_NO_STDC_COMPLEX_KERNELS
  else if constexpr (std::is_same_v<T, std::complex<float>>
                     && has_complex_ufcx_kernels())
  {
    return reinterpret_cast<kernelptr_t<T, U>>(
        integral->tabulate_tensor_complex64);
  }
  else if constexpr (std::is_same_v<T, std::complex<double>>
                     && has_complex_ufcx_kernels())
  {
    return reinterpret_cast<kernelptr_t<T, U>>(
        integral->tabulate_tensor_complex128);
  }
#endif
  else
    throw std::runtime_error("Could not extract kernel from ufcx integral.");
}
 
} // namespace dolfinx::fem::impl