Add device-callable Gauss quadrature to gpu/ headers

GaussLegendre1D: 1-D Gauss-Legendre rules (1-7 points), device-callable.
GaussQuadrature: Full topology dispatcher with n_points/point/weight
methods for EDGE, QUAD, HEX (tensor product), TRI (Dunavant), and
TET (Keast/Walkington). All methods are LIBMESH_DEVICE_INLINE.
fillQuadrature: Host-side convenience wrapper using std::vector.

Author: rochi00

include/Makefile.am

@@ -13,7 +13,9 @@ nobase_include_HEADERS = \
   gpu/kokkos_fe_lagrange_2d.h \
   gpu/kokkos_fe_lagrange_3d.h \
   gpu/kokkos_fe_monomial.h \
-  gpu/kokkos_fe_face_map.h
+  gpu/kokkos_fe_face_map.h \
+  gpu/kokkos_fe_map.h \
+  gpu/kokkos_quadrature.h
 endif
 
 # special handholding for prefix_config.m4 generated files

include/base/libmesh_common.h

@@ -30,6 +30,10 @@
 // The library configuration options
 #include "libmesh/libmesh_config.h"
 
+// Device compilation support — must be included before assert macros
+// so that LIBMESH_DEVICE_ASSERT is available for the Kokkos path.
+#include "libmesh/libmesh_device.h"
+
 // Use actual timestamps or constant dummies (to aid ccache)
 #ifdef LIBMESH_ENABLE_TIMESTAMPS
 #  define  LIBMESH_TIME __TIME__
@@ -287,7 +291,13 @@ extern bool warned_about_auto_ptr;
 #endif
 
 // The libmesh_assert() macro acts like C's assert(), but throws a
-// libmesh_error() (including stack trace, etc) instead of just exiting
+// libmesh_error() (including stack trace, etc) instead of just exiting.
+//
+// In .K translation units (LIBMESH_KOKKOS_COMPILATION defined),
+// LIBMESH_DEVICE_ASSERT is provided by libmesh_device.h using
+// printf + Kokkos::abort() — device-safe across CUDA/HIP/SYCL.
+// The assert macros delegate to it so that both host and device
+// code in the same file get assertion checking.
 #ifdef NDEBUG
 
 #define libmesh_assert_msg(asserted, msg)  ((void) 0)
@@ -299,6 +309,18 @@ extern bool warned_about_auto_ptr;
 #define libmesh_assert_less_equal_msg(expr1,expr2, msg)  ((void) 0)
 #define libmesh_assert_greater_equal_msg(expr1,expr2, msg)  ((void) 0)
 
+#elif defined(LIBMESH_DEVICE_ASSERT)
+
+// Kokkos compilation: use the device-safe assert from libmesh_device.h.
+#define libmesh_assert_msg(asserted, msg)  LIBMESH_DEVICE_ASSERT(asserted)
+#define libmesh_exceptionless_assert_msg(asserted, msg) LIBMESH_DEVICE_ASSERT(asserted)
+#define libmesh_assert_equal_to_msg(expr1,expr2, msg)  LIBMESH_DEVICE_ASSERT((expr1) == (expr2))
+#define libmesh_assert_not_equal_to_msg(expr1,expr2, msg)  LIBMESH_DEVICE_ASSERT((expr1) != (expr2))
+#define libmesh_assert_less_msg(expr1,expr2, msg)  LIBMESH_DEVICE_ASSERT((expr1) < (expr2))
+#define libmesh_assert_greater_msg(expr1,expr2, msg)  LIBMESH_DEVICE_ASSERT((expr1) > (expr2))
+#define libmesh_assert_less_equal_msg(expr1,expr2, msg)  LIBMESH_DEVICE_ASSERT((expr1) <= (expr2))
+#define libmesh_assert_greater_equal_msg(expr1,expr2, msg)  LIBMESH_DEVICE_ASSERT((expr1) >= (expr2))
+
 #else
 
 #define libmesh_assertion_types(expr1,expr2)                            \

include/base/libmesh_device.h

@@ -26,7 +26,24 @@
 // translation units it expands to plain `inline`.
 #ifdef LIBMESH_KOKKOS_COMPILATION
 #  include <Kokkos_Macros.hpp>
+#  include <Kokkos_Abort.hpp>
 #  define LIBMESH_DEVICE_INLINE KOKKOS_INLINE_FUNCTION
+
+// Device-safe assert: uses printf (supported on CUDA/HIP) and
+// Kokkos::abort() for backend-portable device termination.
+// Defined here (not in libmesh_common.h) because Kokkos headers
+// are only available in .K translation units.
+#  ifndef NDEBUG
+#    define LIBMESH_DEVICE_ASSERT(asserted)                             \
+       do { if (!(asserted)) {                                          \
+         printf("libMesh assert failed: %s, file %s, line %d\n",       \
+                #asserted, __FILE__, __LINE__);                         \
+         ::Kokkos::abort("libmesh_assert failed");                     \
+       } } while (0)
+#  else
+#    define LIBMESH_DEVICE_ASSERT(asserted) ((void) 0)
+#  endif
+
 #else
 #  define LIBMESH_DEVICE_INLINE inline
 #endif

include/gpu/kokkos_fe_evaluator.h

@@ -104,6 +104,28 @@ lagrangeTopologyForClassAndOrder(FEElemClass cls, unsigned int order)
 // LAGRANGE_MAP is supported; RATIONAL_BERNSTEIN_MAP requires additional
 // rational-weight data that is not yet threaded through the device path.
 
+// ── Compile-time topology versions (preferred for GPU) ───────────────────
+// Template on FEFamily and ElemType so nvcc only instantiates the specific
+// FEEvaluator specialization. No topology switch means no stack pressure.
+
+/// Compile-time map shape evaluation.
+template <libMesh::FEFamily family, libMesh::ElemType topo>
+LIBMESH_DEVICE_INLINE Real
+nativeMapShape(unsigned int i, Real xi, Real eta, Real zeta)
+{
+  return FEEvaluator<family, topo>::shape(i, xi, eta, zeta);
+}
+
+/// Compile-time map gradient evaluation.
+template <libMesh::FEFamily family, libMesh::ElemType topo>
+LIBMESH_DEVICE_INLINE Real3
+nativeGradMapShape(unsigned int i, Real xi, Real eta, Real zeta)
+{
+  return FEEvaluator<family, topo>::grad_shape(i, xi, eta, zeta);
+}
+
+// ── Runtime topology versions (larger GPU stack usage) ───────────────────
+
 /// Evaluate the i-th geometric map shape function at (xi, eta, zeta).
 LIBMESH_DEVICE_INLINE Real
 nativeMapShape(libMesh::ElemMappingType mapping_type,

include/gpu/kokkos_fe_map.h

@@ -0,0 +1,209 @@
+// Kokkos device-compatible physical map evaluation.
+//
+// All functions are LIBMESH_DEVICE_INLINE — callable from both host and GPU.
+//
+// Two API levels:
+//   1. Template on ElemType (preferred): eliminates the topology switch at
+//      compile time, producing small inlined functions with no stack pressure.
+//   2. Runtime ElemType dispatch: convenient but requires increased CUDA
+//      stack size due to the large switch in nativeMapShape.
+//
+// Given node coordinates and a reference-space point, these functions compute:
+//   - Physical coordinates (xyz)
+//   - Jacobian matrix (reference -> physical)
+//   - Determinant and JxW
+//   - Outward normal (for face integrals)
+
+#pragma once
+
+#include "gpu/kokkos_fe_evaluator.h"
+#include "gpu/kokkos_scalar_types.h"
+
+namespace libMesh::Kokkos
+{
+
+// =========================================================================
+//  Compile-time dispatch (preferred for GPU — no switch overhead)
+//
+//  Template on FEFamily and ElemType so nvcc only instantiates the specific
+//  FEEvaluator specialization. No topology switch means no stack pressure.
+// =========================================================================
+
+/// Compute physical coordinate: xyz = sum_i( phi_i(ref) * node_i )
+template <libMesh::FEFamily family, libMesh::ElemType topo>
+LIBMESH_DEVICE_INLINE Real3
+physicalPoint(const Real3 * nodes,
+              unsigned int n_nodes,
+              Real xi, Real eta, Real zeta)
+{
+  Real3 xyz(0, 0, 0);
+  for (unsigned int i = 0; i < n_nodes; ++i)
+    xyz += nativeMapShape<family, topo>(i, xi, eta, zeta) * nodes[i];
+  return xyz;
+}
+
+/// Compute Jacobian matrix: J_ij = d(x_i)/d(xi_j)
+template <libMesh::FEFamily family, libMesh::ElemType topo>
+LIBMESH_DEVICE_INLINE Real33
+jacobian(const Real3 * nodes,
+         unsigned int n_nodes,
+         Real xi, Real eta, Real zeta)
+{
+  Real33 J(0.0);
+  for (unsigned int k = 0; k < n_nodes; ++k)
+    J += nativeGradMapShape<family, topo>(k, xi, eta, zeta)
+             .cartesian_product(nodes[k]);
+  return J;
+}
+
+/// Compute physical point and Jacobian together (single node loop).
+template <libMesh::FEFamily family, libMesh::ElemType topo>
+LIBMESH_DEVICE_INLINE void
+physicalPointAndJacobian(const Real3 * nodes,
+                         unsigned int n_nodes,
+                         Real xi, Real eta, Real zeta,
+                         Real3 & xyz,
+                         Real33 & J)
+{
+  xyz = Real3(0, 0, 0);
+  J = Real33(0.0);
+  for (unsigned int k = 0; k < n_nodes; ++k)
+  {
+    const Real phi = nativeMapShape<family, topo>(k, xi, eta, zeta);
+    const Real3 grad = nativeGradMapShape<family, topo>(k, xi, eta, zeta);
+    xyz += phi * nodes[k];
+    J += grad.cartesian_product(nodes[k]);
+  }
+}
+
+/// Face Jacobian from side node positions.
+template <libMesh::FEFamily family, libMesh::ElemType face_topo>
+LIBMESH_DEVICE_INLINE Real33
+faceJacobian(const Real3 * face_nodes,
+             unsigned int n_face_nodes,
+             Real xi, Real eta, Real zeta)
+{
+  Real33 J(0.0);
+  for (unsigned int k = 0; k < n_face_nodes; ++k)
+    J += nativeGradMapShape<family, face_topo>(k, xi, eta, zeta)
+             .cartesian_product(face_nodes[k]);
+  return J;
+}
+
+// =========================================================================
+//  Runtime topology dispatch (convenient, but larger GPU stack usage)
+// =========================================================================
+
+/// Compute physical coordinate (runtime topology).
+LIBMESH_DEVICE_INLINE Real3
+physicalPoint(libMesh::ElemMappingType mapping_type,
+              libMesh::ElemType topo,
+              const Real3 * nodes,
+              unsigned int n_nodes,
+              Real xi, Real eta, Real zeta)
+{
+  Real3 xyz(0, 0, 0);
+  for (unsigned int i = 0; i < n_nodes; ++i)
+    xyz += nativeMapShape(mapping_type, topo, i, xi, eta, zeta) * nodes[i];
+  return xyz;
+}
+
+/// Compute Jacobian matrix (runtime topology).
+LIBMESH_DEVICE_INLINE Real33
+jacobian(libMesh::ElemMappingType mapping_type,
+         libMesh::ElemType topo,
+         const Real3 * nodes,
+         unsigned int n_nodes,
+         Real xi, Real eta, Real zeta)
+{
+  Real33 J(0.0);
+  for (unsigned int k = 0; k < n_nodes; ++k)
+    J += nativeGradMapShape(mapping_type, topo, k, xi, eta, zeta)
+             .cartesian_product(nodes[k]);
+  return J;
+}
+
+/// Compute physical point and Jacobian together (runtime topology).
+LIBMESH_DEVICE_INLINE void
+physicalPointAndJacobian(libMesh::ElemMappingType mapping_type,
+                         libMesh::ElemType topo,
+                         const Real3 * nodes,
+                         unsigned int n_nodes,
+                         Real xi, Real eta, Real zeta,
+                         Real3 & xyz,
+                         Real33 & J)
+{
+  xyz = Real3(0, 0, 0);
+  J = Real33(0.0);
+  for (unsigned int k = 0; k < n_nodes; ++k)
+  {
+    const Real phi = nativeMapShape(mapping_type, topo, k, xi, eta, zeta);
+    const Real3 grad = nativeGradMapShape(mapping_type, topo, k, xi, eta, zeta);
+    xyz += phi * nodes[k];
+    J += grad.cartesian_product(nodes[k]);
+  }
+}
+
+/// Face Jacobian (runtime topology).
+LIBMESH_DEVICE_INLINE Real33
+faceJacobian(libMesh::ElemMappingType mapping_type,
+             libMesh::ElemType face_topo,
+             const Real3 * face_nodes,
+             unsigned int n_face_nodes,
+             Real xi, Real eta, Real zeta)
+{
+  Real33 J(0.0);
+  for (unsigned int k = 0; k < n_face_nodes; ++k)
+    J += nativeGradMapShape(mapping_type, face_topo, k, xi, eta, zeta)
+             .cartesian_product(face_nodes[k]);
+  return J;
+}
+
+// =========================================================================
+//  Geometry helpers (topology-independent)
+// =========================================================================
+
+/// Volume element JxW: |det(J)| * quadrature_weight
+LIBMESH_DEVICE_INLINE Real
+volumeJxW(const Real33 & J, unsigned int dim, Real quad_weight)
+{
+  return J.determinant(dim) * quad_weight;
+}
+
+/// Face JxW: surface measure * quadrature_weight
+///   3D: ||J_row0 x J_row1|| * weight
+///   2D: ||J_row0||           * weight
+///   1D: weight (face is a point)
+LIBMESH_DEVICE_INLINE Real
+faceJxW(const Real33 & J, unsigned int parent_dim, Real quad_weight)
+{
+  if (parent_dim == 3)
+    return J.row(0).cross(J.row(1)).norm() * quad_weight;
+  else if (parent_dim == 2)
+    return J.row(0).norm() * quad_weight;
+  else
+    return quad_weight;
+}
+
+/// Outward unit normal from the face Jacobian.
+LIBMESH_DEVICE_INLINE Real3
+faceNormal(const Real33 & J, unsigned int parent_dim)
+{
+  Real3 n(0, 0, 0);
+  if (parent_dim == 3)
+    n = J.row(0).cross(J.row(1));
+  else if (parent_dim == 2)
+  {
+    const Real3 t = J.row(0);
+    n = Real3(-t(1), t(0), 0.0);
+  }
+  else
+    return Real3(1.0, 0.0, 0.0);
+
+  const Real len = n.norm();
+  if (len > 0.0)
+    n *= 1.0 / len;
+  return n;
+}
+
+} // namespace libMesh::Kokkos