From 5a4f31a355907ee1a3058ec00e88dfe24ea0f03c Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Tue, 24 Feb 2026 15:09:28 +0000 Subject: [PATCH 01/17] Implemented MeshManager level versions of get_surface_vertices() + get_surface_connectivity() --- src/mesh_manager_interface.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/mesh_manager_interface.cpp b/src/mesh_manager_interface.cpp index 82cba68b..44c4f6ce 100644 --- a/src/mesh_manager_interface.cpp +++ b/src/mesh_manager_interface.cpp @@ -1,6 +1,8 @@ #include "xdg/mesh_manager_interface.h" #include +#include +#include #include "xdg/config.h" #include "xdg/error.h" From 1c39076a41cc622215e396840cfa54cc058980d6 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Tue, 24 Feb 2026 15:31:14 +0000 Subject: [PATCH 02/17] Refactored get_volume_vertices() and get_volume_connectivity() to be abstracted to MeshManager level --- src/mesh_manager_interface.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mesh_manager_interface.cpp b/src/mesh_manager_interface.cpp index 44c4f6ce..e8d8f009 100644 --- a/src/mesh_manager_interface.cpp +++ b/src/mesh_manager_interface.cpp @@ -301,7 +301,7 @@ std::vector MeshManager::get_surface_vertices(MeshID surface) const return surface_local_mesh_data(surface).vertices; } -std::vector MeshManager::get_surface_connectivity(MeshID surface) const +std::vector MeshManager::get_volume_connectivity(MeshID volume) const { return surface_local_mesh_data(surface).connectivity; } From a17d5db5937990627dd9930b717ebf16b27568ed Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Wed, 25 Feb 2026 18:30:35 +0000 Subject: [PATCH 03/17] Some more cleanup --- src/mesh_manager_interface.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/mesh_manager_interface.cpp b/src/mesh_manager_interface.cpp index e8d8f009..1f192569 100644 --- a/src/mesh_manager_interface.cpp +++ b/src/mesh_manager_interface.cpp @@ -1,8 +1,6 @@ #include "xdg/mesh_manager_interface.h" #include -#include -#include #include "xdg/config.h" #include "xdg/error.h" From c89e4988e24fe8bd6db060c8ff02fcc5ca74dbac Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Fri, 6 Feb 2026 15:33:48 +0000 Subject: [PATCH 04/17] Refactored tests for MOABMeshManager::get_surface_connectivity() and get_volume_connectivity() --- tests/test_moab.cpp | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/tests/test_moab.cpp b/tests/test_moab.cpp index d4abca0c..61f156bc 100644 --- a/tests/test_moab.cpp +++ b/tests/test_moab.cpp @@ -19,6 +19,16 @@ using namespace xdg; using namespace xdg::test; +static bool contains_vertex(const std::vector& vertices, const Vertex& v, double tol) +{ + return std::any_of(vertices.begin(), vertices.end(), + [&](const Vertex& candidate) { + return std::fabs(candidate.x - v.x) <= tol && + std::fabs(candidate.y - v.y) <= tol && + std::fabs(candidate.z - v.z) <= tol; + }); +} + TEST_CASE("Test MOAB Initialization") { std::unique_ptr mesh_manager = std::make_unique(); From af11f4a0837f3f4901a0de3326b90e5355372e47 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Tue, 10 Feb 2026 16:29:11 +0000 Subject: [PATCH 05/17] Added method to return dense array of a volume's connectivity indices --- tests/test_moab.cpp | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/tests/test_moab.cpp b/tests/test_moab.cpp index 61f156bc..d4abca0c 100644 --- a/tests/test_moab.cpp +++ b/tests/test_moab.cpp @@ -19,16 +19,6 @@ using namespace xdg; using namespace xdg::test; -static bool contains_vertex(const std::vector& vertices, const Vertex& v, double tol) -{ - return std::any_of(vertices.begin(), vertices.end(), - [&](const Vertex& candidate) { - return std::fabs(candidate.x - v.x) <= tol && - std::fabs(candidate.y - v.y) <= tol && - std::fabs(candidate.z - v.z) <= tol; - }); -} - TEST_CASE("Test MOAB Initialization") { std::unique_ptr mesh_manager = std::make_unique(); From 4d080ad54e2b0e7496226b388e2bc97902c5f43d Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Wed, 11 Feb 2026 15:12:53 +0000 Subject: [PATCH 06/17] Started on implementing volumetric element tree BVH construction in GPRT --- include/xdg/gprt/ray_tracer.h | 8 +++-- include/xdg/gprt/shared_structs.h | 11 +++++++ src/gprt/dbl_deviceCode.slang | 24 ++++++++++++-- src/gprt/ray_tracer.cpp | 52 ++++++++++++++++++++++++++++--- 4 files changed, 86 insertions(+), 9 deletions(-) diff --git a/include/xdg/gprt/ray_tracer.h b/include/xdg/gprt/ray_tracer.h index 8d24d107..bdeca114 100644 --- a/include/xdg/gprt/ray_tracer.h +++ b/include/xdg/gprt/ray_tracer.h @@ -115,8 +115,9 @@ class GPRTRayTracer : public RayTracer { std::map> rayGenPrograms_; GPRTMissOf missProgram_; - GPRTComputeOf aabbPopulationProgram_; // aabbTriPopulationProgram_; // aabbTetPopulationProgram_; // excludePrimitivesBuffer_; // globalBlasInstances_; // trianglesGeomType_; // tetrahedraGeomType_; // surface_volume_tree_to_accel_map; // Map from XDG::TreeID to GPRTAccel for volume TLAS + std::unordered_map element_volume_tree_to_accel_map; // Map from XDG::TreeID to GPRTAccel for element TLAS + std::vector blas_handles_; // Store BLAS handles so that they can be explicitly referenced in destructor // Global Tree IDs diff --git a/include/xdg/gprt/shared_structs.h b/include/xdg/gprt/shared_structs.h index caa9d553..e1ddf0d0 100644 --- a/include/xdg/gprt/shared_structs.h +++ b/include/xdg/gprt/shared_structs.h @@ -47,6 +47,17 @@ struct DPTriangleGeomData { int num_faces; // Number of faces in the geometry }; +struct DPTetrahedronGeomData { + double3 *vertex; // vertex buffer + float3 *aabbs; // AABB buffer + uint4 *index; // index buffer + int32_t vol_id; + dblRay *ray; // double precision rays + xdg::HitOrientation hitOrientation; + GPRTPrimitiveRef* primitive_refs; + int num_tets; // Number of tetrahedra in the geometry +}; + struct dblRayGenData { dblRay *ray; dblHit *hit; diff --git a/src/gprt/dbl_deviceCode.slang b/src/gprt/dbl_deviceCode.slang index 4e643fd9..55c6bb88 100644 --- a/src/gprt/dbl_deviceCode.slang +++ b/src/gprt/dbl_deviceCode.slang @@ -111,11 +111,11 @@ void point_in_volume(uniform dblRayGenData record, uniform DPTriangleGeomData me } // ------------------------------------------------- Compute Shaders ------------------------------------------------- -/* A shader to compute and store AABB min/maxes in single precision using double precision coords*/ +/* A shader to compute and store AABB min/maxes in single precision using double precision coords for Triangles*/ [shader("compute")] [numthreads(1, 1, 1)] void -populate_aabbs(uint3 DispatchThreadID: SV_DispatchThreadID, uniform DPTriangleGeomData record) { +populate_tri_aabbs(uint3 DispatchThreadID: SV_DispatchThreadID, uniform DPTriangleGeomData record) { int primID = DispatchThreadID.x; int3 indices = record.index[primID]; dp::vec3 A = record.vertex[indices[0]]; @@ -130,6 +130,26 @@ populate_aabbs(uint3 DispatchThreadID: SV_DispatchThreadID, uniform DPTriangleGe record.aabbs[2 * primID + 1] = fpaabbmax; } +/* A shader to compute and store AABB min/maxes in single precision using double precision coords for Tetrahedra*/ +[shader("compute")] +[numthreads(1, 1, 1)] +void +populate_tet_aabbs(uint3 DispatchThreadID: SV_DispatchThreadID, uniform DPTetrahedronGeomData record) { + int primID = DispatchThreadID.x; + int4 indices = record.index[primID]; + double3 A = record.vertex[indices[0]]; + double3 B = record.vertex[indices[1]]; + double3 C = record.vertex[indices[2]]; + double3 D = record.vertex[indices[3]]; + double3 dpaabbmin = min(min(A, B), min(C, D)); + double3 dpaabbmax = max(max(A, B), max(C, D)); + float3 fpaabbmin = float3(dpaabbmin - float3(FLT_EPSILON, FLT_EPSILON, FLT_EPSILON)); + float3 fpaabbmax = float3(dpaabbmax + float3(FLT_EPSILON, FLT_EPSILON, FLT_EPSILON)); + + record.aabbs[2 * primID] = fpaabbmin; + record.aabbs[2 * primID + 1] = fpaabbmax; +} + // ------------------------------------------------ CUSTOM INTERSECTION SHADERS ------------------------------------------------ diff --git a/src/gprt/ray_tracer.cpp b/src/gprt/ray_tracer.cpp index 5ad08e72..d771260a 100644 --- a/src/gprt/ray_tracer.cpp +++ b/src/gprt/ray_tracer.cpp @@ -75,7 +75,8 @@ void GPRTRayTracer::setup_shaders() // TODO: Add Occluded and closest raygen entry points missProgram_ = gprtMissCreate(context_, module_, "ray_fire_miss"); - aabbPopulationProgram_ = gprtComputeCreate(context_, module_, "populate_aabbs"); + aabbTriPopulationProgram_ = gprtComputeCreate(context_, module_, "populate_tri_aabbs"); + aabbTetPopulationProgram_ = gprtComputeCreate(context_, module_, "populate_tet_aabbs"); // Create a "triangle" geometry type and set its closest-hit program trianglesGeomType_ = gprtGeomTypeCreate(context_, GPRT_AABBS); @@ -121,7 +122,7 @@ GPRTRayTracer::create_surface_tree(const std::shared_ptr& mesh_mana auto triangleGeom = gprtGeomCreate(context_, trianglesGeomType_); geom_data = gprtGeomGetParameters(triangleGeom); // pointer to assign data to - // Get storage for vertices + // Get storage for vertices and indices auto vertices = mesh_manager->get_surface_vertices(surf); auto indices = mesh_manager->get_surface_connectivity(surf); std::vector dbl3Vertices; @@ -166,7 +167,7 @@ GPRTRayTracer::create_surface_tree(const std::shared_ptr& mesh_mana geom_data->primitive_refs = gprtBufferGetDevicePointer(primitive_refs_buffer); geom_data->num_faces = num_faces; - gprtComputeLaunch(aabbPopulationProgram_, {num_faces, 1, 1}, {1, 1, 1}, *geom_data); + gprtComputeLaunch(aabbTriPopulationProgram_, {num_faces, 1, 1}, {1, 1, 1}, *geom_data); GPRTAccel blas = gprtAABBAccelCreate(context_, triangleGeom, buildParams_.buildMode); @@ -210,8 +211,49 @@ GPRTRayTracer::create_surface_tree(const std::shared_ptr& mesh_mana ElementTreeID GPRTRayTracer::create_element_tree(const std::shared_ptr& mesh_manager, MeshID volume_id) { - warning("Element trees not currently supported with GPRT ray tracer"); - return TREE_NONE; + auto volume_elements = mesh_manager->get_volume_elements(volume_id); + if (volume_elements.empty()) return TREE_NONE; // No elements in this volume, so no tree to create + + ElementTreeID tree = next_element_tree_id(); + element_trees_.push_back(tree); + + DPTetrahedronGeomData* geom_data = nullptr; + auto tetrahedraGeom = gprtGeomCreate(context_, tetrahedraGeomType_); + + auto vertices = mesh_manager->get_volume_vertices(volume_id); + auto indices = mesh_manager->get_volume_connectivity(volume_id); + std::vector dbl3Vertices; + dbl3Vertices.reserve(vertices.size()); + for (const auto &vertex : vertices) { + dbl3Vertices.push_back({vertex.x, vertex.y, vertex.z}); + } + + // Get storage for indices + std::vector ui4Indices; + ui4Indices.reserve(indices.size() / 4); + for (size_t i = 0; i + 3 < indices.size(); i += 4) { + ui4Indices.emplace_back(indices[i], indices[i + 1], indices[i + 2], indices[i + 3]); + } + + auto vertex_buffer = gprtDeviceBufferCreate(context_, dbl3Vertices.size(), dbl3Vertices.data()); + auto connectivity_buffer = gprtDeviceBufferCreate(context_, ui4Indices.size(), ui4Indices.data()); + auto aabb_buffer = gprtDeviceBufferCreate(context_, 2*volume_elements.size(), 0); // AABBs for each tetrahedron + gprtAABBsSetPositions(tetrahedraGeom, aabb_buffer, volume_elements.size(), 2*sizeof(float3), 0); + + geom_data->aabbs = gprtBufferGetDevicePointer(aabb_buffer); + geom_data->vertex = gprtBufferGetDevicePointer(vertex_buffer); + geom_data->index = gprtBufferGetDevicePointer(connectivity_buffer); + geom_data->num_tets = volume_elements.size(); + geom_data->vol_id = volume_id; + + gprtComputeLaunch(aabbTetPopulationProgram_, {volume_elements.size(), 1, 1}, {1, 1, 1}, *geom_data); + + // GPRTAccel volume_element_accel = nullptr; + + // gprtSolidAccelCreate(context_, volume_element_accel); + + // element_volume_tree_to_accel_map[tree] = volume_element_accel; + }; bool GPRTRayTracer::point_in_volume(SurfaceTreeID tree, From 3eeaab4d4c4835b3e8ea9276e7a6a109e2481e00 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Wed, 11 Feb 2026 16:56:37 +0000 Subject: [PATCH 07/17] Restructured slang shaders in preparation for new tet shaders --- CMakeLists.txt | 13 +- include/xdg/gprt/rt_common.slangh | 53 +++++ include/xdg/shared_enums.h | 1 + src/gprt/dbl_deviceCode.slang | 289 +------------------------- src/gprt/ray_tracer.cpp | 5 +- src/gprt/tetrahedron_rt_shaders.slang | 35 ++++ src/gprt/triangle_rt_shaders.slang | 192 +++++++++++++++++ 7 files changed, 298 insertions(+), 290 deletions(-) create mode 100644 include/xdg/gprt/rt_common.slangh create mode 100644 src/gprt/tetrahedron_rt_shaders.slang create mode 100644 src/gprt/triangle_rt_shaders.slang diff --git a/CMakeLists.txt b/CMakeLists.txt index 8d129141..a9f1afc1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -149,7 +149,14 @@ src/gprt/ray_tracer.cpp list(APPEND xdg_device_codes dbl_deviceCode ) - +list(APPEND xdg_slang_headers +src/gprt/triangle_rt_shaders.slang +src/gprt/tetrahedron_rt_shaders.slang +include/xdg/gprt/shared_structs.h +include/xdg/gprt/rt_common.slangh +include/xdg/shared_enums.h +include/xdg/geometry/dp_math.h +) endif() if (XDG_ENABLE_LIBMESH) @@ -267,9 +274,7 @@ if (XDG_ENABLE_GPRT) OUTPUT_TARGET ${device_code} HEADERS - ${CMAKE_CURRENT_SOURCE_DIR}/include/xdg/gprt/shared_structs.h - ${CMAKE_CURRENT_SOURCE_DIR}/include/xdg/shared_enums.h - ${CMAKE_CURRENT_SOURCE_DIR}/include/xdg/geometry/dp_math.h + ${CMAKE_CURRENT_SOURCE_DIR}/${xdg_slang_headers} SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/src/gprt/${device_code}.slang ) diff --git a/include/xdg/gprt/rt_common.slangh b/include/xdg/gprt/rt_common.slangh new file mode 100644 index 00000000..855b6e3a --- /dev/null +++ b/include/xdg/gprt/rt_common.slangh @@ -0,0 +1,53 @@ +#ifndef XDG_GPRT_RT_COMMON_SLANGH +#define XDG_GPRT_RT_COMMON_SLANGH + +#include "shared_structs.h" +#include "../geometry/dp_math.h" // dp math shared between C++ and slang +#include "../geometry/plucker.h" // Plücker ray-edge intersection test + +// Shared types for both Triangle and Tetrahedron ray tracing pipelines. +[[vk::push_constant]] +dblRayFirePushConstants PC; + +struct DPAttribute +{ + double f64t; // double precision hit distance + int global_prim_id; +}; + +// ------------------------------------------------- Helper functions ------------------------------------------------- + +bool orientation_cull(in double3 ray, in double3 normal, in xdg::HitOrientation orientation) { + if (orientation == xdg::HitOrientation::ANY) return false; // No culling + if (orientation == xdg::HitOrientation::EXITING) return dot(ray, normal) < 0.0; // Cull exiting rays + if (orientation == xdg::HitOrientation::ENTERING) return dot(ray, normal) > 0.0; // Cull entering rays + return false; // Default case, no culling +} + + +/* Function to return the vertex with the lowest coordinates. To force the same + ray-edge computation, the Plücker test needs to use consistent edge + representation. This would be more simple with MOAB handles instead of + coordinates... */ +inline bool first(in double3 a, in double3 b) +{ + if (a[0] < b[0]) return true; + + if (a[0] == b[0] && a[1] < b[1]) return true; + + if (a[1] == b[1] && a[2] < b[2]) return true; + + return false; +} + +float next_after(float a) { + uint a_ = asuint(a); + if (a < 0) { + a_--; + } else { + a_++; + } + return asfloat(a_); +} + +#endif // XDG_GPRT_RT_COMMON_SLANGH diff --git a/include/xdg/shared_enums.h b/include/xdg/shared_enums.h index f6198f60..a20c27fa 100644 --- a/include/xdg/shared_enums.h +++ b/include/xdg/shared_enums.h @@ -4,6 +4,7 @@ namespace xdg { enum PointInVolume : int { + UNSET = -1, OUTSIDE = 0, INSIDE = 1 }; diff --git a/src/gprt/dbl_deviceCode.slang b/src/gprt/dbl_deviceCode.slang index 55c6bb88..49323e98 100644 --- a/src/gprt/dbl_deviceCode.slang +++ b/src/gprt/dbl_deviceCode.slang @@ -1,286 +1,7 @@ -#include "../../include/xdg/gprt/shared_structs.h" -#include "../../include/xdg/geometry/plucker.h" +#ifndef DBL_DEVICE_CODE_SLANGH +#define DBL_DEVICE_CODE_SLANGH -/* -For now we have to use relative paths for includes, which is not ideal. If https://github.com/gprt-org/GPRT/pull/82 gets -merged into GPRT we will be able to more robustly include these headers in the manner below: +#include "triangle_rt_shaders.slang" +#include "tetrahedron_rt_shaders.slang" -#include "xdg/gprt/shared_structs.h" -#include "xdg/geometry/plucker.h" -*/ - -[[vk::push_constant]] -dblRayFirePushConstants PC; - -struct RayFirePayload { - double distance; // Distance to intersection - int surf_id; // ID of the surface hit - SurfaceAccelerationStructure tlas; - int primitive_id; // ID of the primitive hit - xdg::PointInVolume piv; // Point in volume check (0 for outside, 1 for inside) -}; - - -struct DPAttribute -{ - double f64t; // double precision hit distance - int global_prim_id; -}; - - -[shader("closesthit")] -void ray_fire_hit(uniform DPTriangleGeomData record, inout RayFirePayload payload, in DPAttribute attr) { - // Distance from the ray origin to the hit point - uint hit_kind = HitKind(); - uint rayID = DispatchRaysIndex().x; - - // There is some logic for handling next volumes inside the h5m-reader which I could make use of too - // TODO : Should the dblHit struct return the next volume ID for the ray back to the host - - payload.piv = (hit_kind == HIT_KIND_TRIANGLE_FRONT_FACE) - ? xdg::PointInVolume::OUTSIDE - : xdg::PointInVolume::INSIDE; - - int instanceID = InstanceID(); - - payload.distance = attr.f64t; - payload.surf_id = record.surf_id; - payload.primitive_id = attr.global_prim_id; -} - -[shader("miss")] -void ray_fire_miss(inout RayFirePayload payload) { - // Set the miss payload to default values - payload.distance = -1.0f; - payload.surf_id = -1; - payload.primitive_id = -1; -} - -// This ray generation program will kick off the ray tracing process, -// generating rays and tracing them into the world. -[shader("raygeneration")] -void ray_fire(uniform dblRayGenData record, uniform DPTriangleGeomData mesh) { - RayFirePayload payload; - uint rayID = DispatchRaysIndex().x; - - // Trace the ray into the scene - RayDesc rayDesc; - rayDesc.Origin = float3(record.ray[rayID].origin); - rayDesc.Direction = normalize(float3(record.ray[rayID].direction)); - rayDesc.TMin = float(record.ray[rayID].tMin); - rayDesc.TMax = float(record.ray[rayID].tMax); - - SurfaceAccelerationStructure world = record.ray[rayID].volume_accel; - - // Pass the ray's origin and direction to the payload - payload.distance = -1.0f; - payload.surf_id = -1; - payload.tlas = world; - - TraceRay(world, RAY_FLAG_NONE, 0xff, 0, 1, rayDesc, payload); - - // Store the distance to the hit point and the surface ID in buffers for CPU - record.hit[rayID].distance = payload.distance; - record.hit[rayID].surf_id = payload.surf_id; - record.hit[rayID].primitive_id = payload.primitive_id; -} - -[shader("raygeneration")] -void point_in_volume(uniform dblRayGenData record, uniform DPTriangleGeomData mesh) { - RayFirePayload payload; - uint rayID = DispatchRaysIndex().x; - - // Trace the ray into the scene - RayDesc rayDesc; - rayDesc.Origin = float3(record.ray[rayID].origin); - rayDesc.Direction = float3(normalize(record.ray[rayID].direction)); - rayDesc.TMin = float(record.ray[rayID].tMin); - rayDesc.TMax = float(record.ray[rayID].tMax); - - SurfaceAccelerationStructure world = record.ray[rayID].volume_accel; - - // Pass the ray's origin and direction to the payload - payload.surf_id = -1; - payload.tlas = world; - payload.piv = xdg::PointInVolume::OUTSIDE; // Initialize point in volume check result to outside (0) - - TraceRay(world, RAY_FLAG_NONE, 0xff, 0, 1, rayDesc, payload); - - record.hit[rayID].surf_id = payload.surf_id; - record.hit[rayID].piv = payload.piv; // Point in volume check result -} - -// ------------------------------------------------- Compute Shaders ------------------------------------------------- -/* A shader to compute and store AABB min/maxes in single precision using double precision coords for Triangles*/ -[shader("compute")] -[numthreads(1, 1, 1)] -void -populate_tri_aabbs(uint3 DispatchThreadID: SV_DispatchThreadID, uniform DPTriangleGeomData record) { - int primID = DispatchThreadID.x; - int3 indices = record.index[primID]; - dp::vec3 A = record.vertex[indices[0]]; - dp::vec3 B = record.vertex[indices[1]]; - dp::vec3 C = record.vertex[indices[2]]; - dp::vec3 dpaabbmin = min(min(A, B), C); - dp::vec3 dpaabbmax = max(max(A, B), C); - float3 fpaabbmin = float3(dpaabbmin - float3(FLT_EPSILON, FLT_EPSILON, FLT_EPSILON)); - float3 fpaabbmax = float3(dpaabbmax + float3(FLT_EPSILON, FLT_EPSILON, FLT_EPSILON)); - - record.aabbs[2 * primID] = fpaabbmin; - record.aabbs[2 * primID + 1] = fpaabbmax; -} - -/* A shader to compute and store AABB min/maxes in single precision using double precision coords for Tetrahedra*/ -[shader("compute")] -[numthreads(1, 1, 1)] -void -populate_tet_aabbs(uint3 DispatchThreadID: SV_DispatchThreadID, uniform DPTetrahedronGeomData record) { - int primID = DispatchThreadID.x; - int4 indices = record.index[primID]; - double3 A = record.vertex[indices[0]]; - double3 B = record.vertex[indices[1]]; - double3 C = record.vertex[indices[2]]; - double3 D = record.vertex[indices[3]]; - double3 dpaabbmin = min(min(A, B), min(C, D)); - double3 dpaabbmax = max(max(A, B), max(C, D)); - float3 fpaabbmin = float3(dpaabbmin - float3(FLT_EPSILON, FLT_EPSILON, FLT_EPSILON)); - float3 fpaabbmax = float3(dpaabbmax + float3(FLT_EPSILON, FLT_EPSILON, FLT_EPSILON)); - - record.aabbs[2 * primID] = fpaabbmin; - record.aabbs[2 * primID + 1] = fpaabbmax; -} - -// ------------------------------------------------ CUSTOM INTERSECTION SHADERS ------------------------------------------------ - - -/* 1D ray generation intersection with a double precision triangle using the Plucker intersection algorithm*/ -[shader("intersection")] -void DPTrianglePluckerIntersection(uniform DPTriangleGeomData record) -{ - int primID = PrimitiveIndex(); - int global_prim_id = record.primitive_refs[primID].id; - - uint rayID = DispatchRaysIndex().x; - uint nRays = DispatchRaysDimensions().x; - - if (rayID >= nRays) return; - - // Load vertices - int3 indices = record.index[primID]; - dp::vec3 v0 = record.vertex[indices[0]]; - dp::vec3 v1 = record.vertex[indices[1]]; - dp::vec3 v2 = record.vertex[indices[2]]; - dp::vec3 vertices[3] = { v0, v1, v2 }; - - // Load rays - dp::vec3 origin = record.ray[rayID].origin; - dp::vec3 direction = record.ray[rayID].direction; - double tMin = record.ray[rayID].tMin; - double tMax = record.ray[rayID].tMax; - - bool useOrientation = false; - int orientation = 0; - - xdg::PluckerIntersectionResult result = xdg::plucker_ray_tri_intersect(vertices, - origin, - direction, - tMax, - tMin, - useOrientation, - orientation); - if (result.hit == false) return; // No intersection - double t = result.t; - - DPAttribute attr; - attr.f64t = t; - - float f32t = float(t); - if (double(f32t) < t) f32t = next_after(f32t); - - dp::vec3 norm = record.normals[primID]; // recover double precision normal. TODO - Should we calculate from vertices instead? - - // sense adjustment of normal - if (record.ray[rayID].volume_tree == record.reverse_tree) - { - norm = -norm; - } - - double norm_dot_dir = dp::dot(norm, direction); - uint hit_kind = norm_dot_dir < 0 ? HIT_KIND_TRIANGLE_FRONT_FACE - : HIT_KIND_TRIANGLE_BACK_FACE; - - xdg::HitOrientation hitOrientation = record.ray[rayID].hitOrientation; - - if (orientation_cull(direction, norm, hitOrientation)) - { - return; - } - - for (int i = 0; i < record.ray[rayID].exclude_count; ++i) - { - if (record.ray[rayID].exclude_primitives[i] == global_prim_id) { - return; - } - } - attr.global_prim_id = global_prim_id; - ReportHit(f32t, hit_kind, attr); -} - - -// ------------------------------------------------- Helper functions ------------------------------------------------- - -bool orientation_cull(in dp::vec3 ray, in dp::vec3 normal, in xdg::HitOrientation orientation) { - if (orientation == xdg::HitOrientation::ANY) return false; // No culling - if (orientation == xdg::HitOrientation::EXITING) return dp::dot(ray, normal) < 0.0; // Cull exiting rays - if (orientation == xdg::HitOrientation::ENTERING) return dp::dot(ray, normal) > 0.0; // Cull entering rays - return false; // Default case, no culling -} - -// Plucker coordinate -double plucker_edge_test(in dp::vec3 vertexa, in dp::vec3 vertexb, in dp::vec3 ray, in dp::vec3 ray_normal) -{ - double pip; - const double near_zero = 10 * DBL_EPSILON; - - if (first(vertexa, vertexb)) - { - dp::vec3 edge = vertexb - vertexa; - dp::vec3 edge_normal = dp::cross(edge, vertexa); - pip = dp::dot(ray, edge_normal) + dp::dot(ray_normal, edge); - } - else - { - dp::vec3 edge = vertexa - vertexb; - dp::vec3 edge_normal = dp::cross(edge, vertexb); - pip = dp::dot(ray, edge_normal) + dp::dot(ray_normal, edge); - pip = -pip; - } - - if (near_zero > abs(pip)) pip = 0.0; - return pip; -} - -/* Function to return the vertex with the lowest coordinates. To force the same - ray-edge computation, the Plücker test needs to use consistent edge - representation. This would be more simple with MOAB handles instead of - coordinates... */ -inline bool first(in dp::vec3 a, in dp::vec3 b) -{ - if (a[0] < b[0]) return true; - - if (a[0] == b[0] && a[1] < b[1]) return true; - - if (a[1] == b[1] && a[2] < b[2]) return true; - - return false; -} - -float next_after(float a) { - uint a_ = asuint(a); - if (a < 0) { - a_--; - } else { - a_++; - } - return asfloat(a_); -} \ No newline at end of file +#endif // DBL_DEVICE_CODE_SLANGH \ No newline at end of file diff --git a/src/gprt/ray_tracer.cpp b/src/gprt/ray_tracer.cpp index d771260a..dc22d5d9 100644 --- a/src/gprt/ray_tracer.cpp +++ b/src/gprt/ray_tracer.cpp @@ -248,12 +248,13 @@ GPRTRayTracer::create_element_tree(const std::shared_ptr& mesh_mana gprtComputeLaunch(aabbTetPopulationProgram_, {volume_elements.size(), 1, 1}, {1, 1, 1}, *geom_data); - // GPRTAccel volume_element_accel = nullptr; + GPRTAccel volume_element_accel = gprtAABBAccelCreate(context_, tetrahedraGeom, buildParams_.buildMode); - // gprtSolidAccelCreate(context_, volume_element_accel); + gprtAccelBuild(context_, volume_element_accel, buildParams_); // element_volume_tree_to_accel_map[tree] = volume_element_accel; + return tree; }; bool GPRTRayTracer::point_in_volume(SurfaceTreeID tree, diff --git a/src/gprt/tetrahedron_rt_shaders.slang b/src/gprt/tetrahedron_rt_shaders.slang new file mode 100644 index 00000000..d862fdc5 --- /dev/null +++ b/src/gprt/tetrahedron_rt_shaders.slang @@ -0,0 +1,35 @@ +#ifndef XDG_GPRT_TETRAHEDRON_RT_SHADERS_SLANG +#define XDG_GPRT_TETRAHEDRON_RT_SHADERS_SLANG + +#include "../../include/xdg/gprt/rt_common.slangh" + +struct SolidRayFirePayload { + double distance; // Distance to intersection + int surf_id; // ID of the surface hit + SolidAccelerationStructure tlas; + int primitive_id; // ID of the primitive hit +}; + + +// ------------------------------------------------- Compute Shaders ------------------------------------------------- +/* A shader to compute and store AABB min/maxes in single precision using double precision coords for Tetrahedra*/ +[shader("compute")] +[numthreads(1, 1, 1)] +void +populate_tet_aabbs(uint3 DispatchThreadID: SV_DispatchThreadID, uniform DPTetrahedronGeomData record) { + int primID = DispatchThreadID.x; + int4 indices = record.index[primID]; + double3 A = record.vertex[indices[0]]; + double3 B = record.vertex[indices[1]]; + double3 C = record.vertex[indices[2]]; + double3 D = record.vertex[indices[3]]; + double3 dpaabbmin = min(min(A, B), min(C, D)); + double3 dpaabbmax = max(max(A, B), max(C, D)); + float3 fpaabbmin = float3(dpaabbmin - float3(FLT_EPSILON, FLT_EPSILON, FLT_EPSILON)); + float3 fpaabbmax = float3(dpaabbmax + float3(FLT_EPSILON, FLT_EPSILON, FLT_EPSILON)); + + record.aabbs[2 * primID] = fpaabbmin; + record.aabbs[2 * primID + 1] = fpaabbmax; +} + +#endif // XDG_GPRT_TETRAHEDRON_RT_SHADERS_SLANG diff --git a/src/gprt/triangle_rt_shaders.slang b/src/gprt/triangle_rt_shaders.slang new file mode 100644 index 00000000..b0d90dfa --- /dev/null +++ b/src/gprt/triangle_rt_shaders.slang @@ -0,0 +1,192 @@ +#ifndef XDG_GPRT_TRIANGLE_RT_SHADERS_SLANG +#define XDG_GPRT_TRIANGLE_RT_SHADERS_SLANG + +#include "../../include/xdg/gprt/rt_common.slangh" + +struct SurfaceRayFirePayload { + double distance; // Distance to intersection + int surf_id; // ID of the surface hit + SurfaceAccelerationStructure tlas; + int primitive_id; // ID of the primitive hit + xdg::PointInVolume piv; // Point in volume check (0 for outside, 1 for inside) +}; + +[shader("closesthit")] +void ray_fire_hit(uniform DPTriangleGeomData record, inout SurfaceRayFirePayload payload, in DPAttribute attr) { + // Distance from the ray origin to the hit point + uint hit_kind = HitKind(); + uint rayID = DispatchRaysIndex().x; + + // There is some logic for handling next volumes inside the h5m-reader which I could make use of too + // TODO : Should the dblHit struct return the next volume ID for the ray back to the host + + payload.piv = (hit_kind == HIT_KIND_TRIANGLE_FRONT_FACE) + ? xdg::PointInVolume::OUTSIDE + : xdg::PointInVolume::INSIDE; + + int instanceID = InstanceID(); + + payload.distance = attr.f64t; + payload.surf_id = record.surf_id; + payload.primitive_id = attr.global_prim_id; +} + +[shader("miss")] +void ray_fire_miss(inout SurfaceRayFirePayload payload) { + // Set the miss payload to default values + payload.distance = -1.0f; + payload.surf_id = -1; + payload.primitive_id = -1; + payload.piv = xdg::PointInVolume::UNSET; +} + +// This ray generation program will kick off the ray tracing process, +// generating rays and tracing them into the world. +[shader("raygeneration")] +void ray_fire(uniform dblRayGenData record, uniform DPTriangleGeomData mesh) { + SurfaceRayFirePayload payload; + uint rayID = DispatchRaysIndex().x; + + // Trace the ray into the scene + RayDesc rayDesc; + rayDesc.Origin = float3(record.ray[rayID].origin); + rayDesc.Direction = normalize(float3(record.ray[rayID].direction)); + rayDesc.TMin = float(record.ray[rayID].tMin); + rayDesc.TMax = float(record.ray[rayID].tMax); + + SurfaceAccelerationStructure world = record.ray[rayID].volume_accel; + + // Pass the ray's origin and direction to the payload + payload.distance = -1.0f; + payload.surf_id = -1; + payload.tlas = world; + + TraceRay(world, RAY_FLAG_NONE, 0xff, 0, 1, rayDesc, payload); + + // Store the distance to the hit point and the surface ID in buffers for CPU + record.hit[rayID].distance = payload.distance; + record.hit[rayID].surf_id = payload.surf_id; + record.hit[rayID].primitive_id = payload.primitive_id; +} + +[shader("raygeneration")] +void point_in_volume(uniform dblRayGenData record, uniform DPTriangleGeomData mesh) { + SurfaceRayFirePayload payload; + uint rayID = DispatchRaysIndex().x; + + // Trace the ray into the scene + RayDesc rayDesc; + rayDesc.Origin = float3(record.ray[rayID].origin); + rayDesc.Direction = float3(normalize(record.ray[rayID].direction)); + rayDesc.TMin = float(record.ray[rayID].tMin); + rayDesc.TMax = float(record.ray[rayID].tMax); + + SurfaceAccelerationStructure world = record.ray[rayID].volume_accel; + + // Pass the ray's origin and direction to the payload + payload.surf_id = -1; + payload.tlas = world; + payload.piv = xdg::PointInVolume::OUTSIDE; // Initialize point in volume check result to outside (0) + + TraceRay(world, RAY_FLAG_NONE, 0xff, 0, 1, rayDesc, payload); + + record.hit[rayID].surf_id = payload.surf_id; + record.hit[rayID].piv = payload.piv; // Point in volume check result +} + +// ------------------------------------------------- Compute Shaders ------------------------------------------------- +/* A shader to compute and store AABB min/maxes in single precision using double precision coords for Triangles*/ +[shader("compute")] +[numthreads(1, 1, 1)] +void +populate_tri_aabbs(uint3 DispatchThreadID: SV_DispatchThreadID, uniform DPTriangleGeomData record) { + int primID = DispatchThreadID.x; + int3 indices = record.index[primID]; + double3 A = record.vertex[indices[0]]; + double3 B = record.vertex[indices[1]]; + double3 C = record.vertex[indices[2]]; + double3 dpaabbmin = min(min(A, B), C); + double3 dpaabbmax = max(max(A, B), C); + float3 fpaabbmin = float3(dpaabbmin - float3(FLT_EPSILON, FLT_EPSILON, FLT_EPSILON)); + float3 fpaabbmax = float3(dpaabbmax + float3(FLT_EPSILON, FLT_EPSILON, FLT_EPSILON)); + + record.aabbs[2 * primID] = fpaabbmin; + record.aabbs[2 * primID + 1] = fpaabbmax; +} + +// ------------------------------------------------ CUSTOM INTERSECTION SHADERS ------------------------------------------------ + +/* 1D ray generation intersection with a double precision triangle using the Plucker intersection algorithm*/ +[shader("intersection")] +void DPTrianglePluckerIntersection(uniform DPTriangleGeomData record) +{ + int primID = PrimitiveIndex(); + int global_prim_id = record.primitive_refs[primID].id; + + uint rayID = DispatchRaysIndex().x; + uint nRays = DispatchRaysDimensions().x; + + if (rayID >= nRays) return; + + // Load vertices + int3 indices = record.index[primID]; + dp::vec3 v0 = record.vertex[indices[0]]; + dp::vec3 v1 = record.vertex[indices[1]]; + dp::vec3 v2 = record.vertex[indices[2]]; + dp::vec3 vertices[3] = { v0, v1, v2 }; + + // Load rays + dp::vec3 origin = record.ray[rayID].origin; + dp::vec3 direction = record.ray[rayID].direction; + double tMin = record.ray[rayID].tMin; + double tMax = record.ray[rayID].tMax; + + bool useOrientation = false; + int orientation = 0; + + xdg::PluckerIntersectionResult result = xdg::plucker_ray_tri_intersect(vertices, + origin, + direction, + tMax, + tMin, + useOrientation, + orientation); + if (result.hit == false) return; // No intersection + double t = result.t; + + DPAttribute attr; + attr.f64t = t; + + float f32t = float(t); + if (double(f32t) < t) f32t = next_after(f32t); + + dp::vec3 norm = record.normals[primID]; // recover double precision normal. TODO - Should we calculate from vertices instead? + + // sense adjustment of normal + if (record.ray[rayID].volume_tree == record.reverse_tree) + { + norm = -norm; + } + + double norm_dot_dir = dp::dot(norm, direction); + uint hit_kind = norm_dot_dir < 0 ? HIT_KIND_TRIANGLE_FRONT_FACE + : HIT_KIND_TRIANGLE_BACK_FACE; + + xdg::HitOrientation hitOrientation = record.ray[rayID].hitOrientation; + + if (orientation_cull(direction, norm, hitOrientation)) + { + return; + } + + for (int i = 0; i < record.ray[rayID].exclude_count; ++i) + { + if (record.ray[rayID].exclude_primitives[i] == global_prim_id) { + return; + } + } + attr.global_prim_id = global_prim_id; + ReportHit(f32t, hit_kind, attr); +} + +#endif // XDG_GPRT_TRIANGLE_RT_SHADERS_SLANG From afb63349bb156ab8d026b1805295ca074a142c05 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Fri, 13 Feb 2026 17:58:34 +0000 Subject: [PATCH 08/17] Added the required shaders for find_element on GPU --- include/xdg/gprt/rt_common.slangh | 66 ++++++++++++++++++++++++ include/xdg/gprt/shared_structs.h | 3 +- src/gprt/ray_tracer.cpp | 4 +- src/gprt/tetrahedron_rt_shaders.slang | 72 +++++++++++++++++++++++++-- src/gprt/triangle_rt_shaders.slang | 4 +- 5 files changed, 140 insertions(+), 9 deletions(-) diff --git a/include/xdg/gprt/rt_common.slangh b/include/xdg/gprt/rt_common.slangh index 855b6e3a..80d7ac70 100644 --- a/include/xdg/gprt/rt_common.slangh +++ b/include/xdg/gprt/rt_common.slangh @@ -15,6 +15,10 @@ struct DPAttribute int global_prim_id; }; +static const double PLUCKER_ZERO_TOL = double(20.0) * DBL_EPSILON; + +double3x3 inverse3x3(in double3x3 matrix); + // ------------------------------------------------- Helper functions ------------------------------------------------- bool orientation_cull(in double3 ray, in double3 normal, in xdg::HitOrientation orientation) { @@ -50,4 +54,66 @@ float next_after(float a) { return asfloat(a_); } +// slang doesnt seem to have an inverse matrix intrinsic so I needed to write my own unfortunately :/ +double3x3 inverse3x3(in double3x3 matrix) +{ + double m00 = matrix[0][0], m01 = matrix[0][1], m02 = matrix[0][2]; + double m10 = matrix[1][0], m11 = matrix[1][1], m12 = matrix[1][2]; + double m20 = matrix[2][0], m21 = matrix[2][1], m22 = matrix[2][2]; + + // Cofactor matrix. + double c00 = (m11 * m22 - m12 * m21); + double c01 = -(m10 * m22 - m12 * m20); + double c02 = (m10 * m21 - m11 * m20); + double c10 = -(m01 * m22 - m02 * m21); + double c11 = (m00 * m22 - m02 * m20); + double c12 = -(m00 * m21 - m01 * m20); + double c20 = (m01 * m12 - m02 * m11); + double c21 = -(m00 * m12 - m02 * m10); + double c22 = (m00 * m11 - m01 * m10); + + double det = m00 * c00 + m01 * c01 + m02 * c02; + + // Adjugate is transpose(cofactor). + double3x3 adj = { + { c00, c10, c20 }, + { c01, c11, c21 }, + { c02, c12, c22 } + }; + return adj * (double(1.0) / det); +} + +bool plucker_tet_containment_test(in double3 point, + in double3 v0, + in double3 v1, + in double3 v2, + in double3 v3) { + // Create matrix T = [v1 - v0, v2 - v0, v3 - v0] + double3 e0 = v1 - v0; + double3 e1 = v2 - v0; + double3 e2 = v3 - v0; + double3x3 T = { { e0.x, e0.y, e0.z }, + { e1.x, e1.y, e1.z }, + { e2.x, e2.y, e2.z } }; + + // Vector from v0 to point + double3 rhs = point - v0; + + // Solve T * [λ1, λ2, λ3] = rhs + double3 lambda123 = mul(inverse3x3(T), double3(rhs.x, rhs.y, rhs.z)); + + // Compute λ0 + double lambda0 = double(1.0) - (lambda123.x + lambda123.y + lambda123.z); + + // Final barycentric coordinate vector + double4 bary = { lambda0, lambda123.x, lambda123.y, lambda123.z }; + + // Check all λ_i in [0, 1] + for (int i = 0; i < 4; ++i) { + if (bary[i] < -PLUCKER_ZERO_TOL || bary[i] > double(1.0) + PLUCKER_ZERO_TOL) + return false; + } + return true; +} + #endif // XDG_GPRT_RT_COMMON_SLANGH diff --git a/include/xdg/gprt/shared_structs.h b/include/xdg/gprt/shared_structs.h index e1ddf0d0..f7aa66fa 100644 --- a/include/xdg/gprt/shared_structs.h +++ b/include/xdg/gprt/shared_structs.h @@ -18,7 +18,8 @@ struct dblRay int32_t exclude_count; // Number of excluded primitives xdg::HitOrientation hitOrientation; int volume_tree; // TreeID of the volume being queried - SurfaceAccelerationStructure volume_accel; // The volume accel + SurfaceAccelerationStructure volume_accel_surf; // The volume accel for surface acceleration structure + SolidAccelerationStructure volume_accel_solid; // The volume accel for solid acceleration structure }; struct dblHit diff --git a/src/gprt/ray_tracer.cpp b/src/gprt/ray_tracer.cpp index dc22d5d9..a4116563 100644 --- a/src/gprt/ray_tracer.cpp +++ b/src/gprt/ray_tracer.cpp @@ -272,7 +272,7 @@ bool GPRTRayTracer::point_in_volume(SurfaceTreeID tree, gprtBufferMap(rayHitBuffers_.ray); // Update the ray input buffer dblRay* ray = gprtBufferGetHostPointer(rayHitBuffers_.ray); - ray[0].volume_accel = gprtAccelGetDeviceAddress(volume); + ray[0].volume_accel_surf = gprtAccelGetDeviceAddress(volume); ray[0].origin = {point.x, point.y, point.z}; ray[0].direction = {directionUsed.x, directionUsed.y, directionUsed.z}; ray[0].tMax = INFTY; // Set a large distance limit @@ -329,7 +329,7 @@ std::pair GPRTRayTracer::ray_fire(SurfaceTreeID tree, gprtBufferMap(rayHitBuffers_.ray); // Update the ray input buffer dblRay* ray = gprtBufferGetHostPointer(rayHitBuffers_.ray); - ray[0].volume_accel = gprtAccelGetDeviceAddress(volume); + ray[0].volume_accel_surf = gprtAccelGetDeviceAddress(volume); ray[0].origin = {origin.x, origin.y, origin.z}; ray[0].direction = {direction.x, direction.y, direction.z}; ray[0].tMax = dist_limit; diff --git a/src/gprt/tetrahedron_rt_shaders.slang b/src/gprt/tetrahedron_rt_shaders.slang index d862fdc5..d9e9cc95 100644 --- a/src/gprt/tetrahedron_rt_shaders.slang +++ b/src/gprt/tetrahedron_rt_shaders.slang @@ -4,10 +4,7 @@ #include "../../include/xdg/gprt/rt_common.slangh" struct SolidRayFirePayload { - double distance; // Distance to intersection - int surf_id; // ID of the surface hit - SolidAccelerationStructure tlas; - int primitive_id; // ID of the primitive hit + int primitiveid; // ID of the primitive the point is contained within }; @@ -32,4 +29,71 @@ populate_tet_aabbs(uint3 DispatchThreadID: SV_DispatchThreadID, uniform DPTetrah record.aabbs[2 * primID + 1] = fpaabbmax; } +[shader("raygeneration")] +void find_element(uniform dblRayGenData record, uniform DPTetrahedronGeomData mesh) { + // This shader is only used to find the element containing a point, so we can ignore the raygen data and just use the ray's origin as the point to test for containment + SolidRayFirePayload payload; + uint rayID = DispatchRaysIndex().x; + PointDesc pointDesc; + pointDesc.Origin = float3(mesh.ray[rayID].origin); + + SolidAccelerationStructure world = mesh.ray[rayID].volume_accel_solid; + + // Pass the ray's origin and direction to the payload + payload.primitiveid = -1; // Initialize to -1 to indicate the point is + + TracePoint(world, RAY_FLAG_NONE, 0xff, 0, 1, pointDesc, payload); + + // Store the ID of the containing element in the hit buffer for CPU readback + record.hit[rayID].primitive_id = payload.primitiveid; +} + +[shader("miss")] +void tet_miss(uniform DPTetrahedronGeomData record, inout SolidRayFirePayload payload) { + payload.primitiveid = -1; // Mark the ray as outside all tetrahedra +} + +[shader("closesthit")] +void tet_contain_hit(uniform DPTetrahedronGeomData record, inout SolidRayFirePayload payload, in DPAttribute attr) { + payload.primitiveid = attr.global_prim_id; // Store the ID of the tet that contains the point +} + + +/* 1D ray generation intersection with a double precision Tetrahedron using the Plucker intersection algorithm*/ +[shader("intersection")] +void DPTetrahedronPluckerIntersection(uniform DPTetrahedronGeomData record) +{ + int primID = PrimitiveIndex(); + int global_prim_id = record.primitive_refs[primID].id; + + uint rayID = DispatchRaysIndex().x; + uint nRays = DispatchRaysDimensions().x; + uint flags = RayFlags(); + + if (rayID >= nRays) { + return; + } + + int4 indices = record.index[primID]; + double3 v0 = record.vertex[indices[0]]; + double3 v1 = record.vertex[indices[1]]; + double3 v2 = record.vertex[indices[2]]; + double3 v3 = record.vertex[indices[3]]; + + double3 origin = record.ray[rayID].origin; + + bool inside = plucker_tet_containment_test(origin, v0, v1, v2, v3); + + if (!inside) return; + + // If we are inside the tet we need to report a hit to move through the rest of the RT pipeline + float f32t = 0.0f; // We dont care about intersection distance + uint hit_kind = 0; // No hit kind since it doesn't matter + + DPAttribute attr; + attr.global_prim_id = global_prim_id; + + ReportHit(f32t, hit_kind, attr); +} + #endif // XDG_GPRT_TETRAHEDRON_RT_SHADERS_SLANG diff --git a/src/gprt/triangle_rt_shaders.slang b/src/gprt/triangle_rt_shaders.slang index b0d90dfa..641f4077 100644 --- a/src/gprt/triangle_rt_shaders.slang +++ b/src/gprt/triangle_rt_shaders.slang @@ -54,7 +54,7 @@ void ray_fire(uniform dblRayGenData record, uniform DPTriangleGeomData mesh) { rayDesc.TMin = float(record.ray[rayID].tMin); rayDesc.TMax = float(record.ray[rayID].tMax); - SurfaceAccelerationStructure world = record.ray[rayID].volume_accel; + SurfaceAccelerationStructure world = record.ray[rayID].volume_accel_surf; // Pass the ray's origin and direction to the payload payload.distance = -1.0f; @@ -81,7 +81,7 @@ void point_in_volume(uniform dblRayGenData record, uniform DPTriangleGeomData me rayDesc.TMin = float(record.ray[rayID].tMin); rayDesc.TMax = float(record.ray[rayID].tMax); - SurfaceAccelerationStructure world = record.ray[rayID].volume_accel; + SurfaceAccelerationStructure world = record.ray[rayID].volume_accel_surf; // Pass the ray's origin and direction to the payload payload.surf_id = -1; From af36622238a58156ecd9b7e5ccbffe3f6eceb17d Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Thu, 19 Feb 2026 14:44:48 +0000 Subject: [PATCH 09/17] Added rest of plumbing to make find_element work with GPRT and 1 ray --- include/xdg/gprt/ray_tracer.h | 21 +++--- include/xdg/gprt/rt_common.slangh | 6 +- include/xdg/gprt/shared_structs.h | 8 +++ src/gprt/ray_tracer.cpp | 92 +++++++++++++++++++++++---- src/gprt/tetrahedron_rt_shaders.slang | 15 ++--- src/gprt/triangle_rt_shaders.slang | 8 +-- 6 files changed, 111 insertions(+), 39 deletions(-) diff --git a/include/xdg/gprt/ray_tracer.h b/include/xdg/gprt/ray_tracer.h index bdeca114..4e759fe6 100644 --- a/include/xdg/gprt/ray_tracer.h +++ b/include/xdg/gprt/ray_tracer.h @@ -22,7 +22,8 @@ enum class RayGenType { RAY_FIRE, POINT_IN_VOLUME, OCCLUDED, - CLOSEST + CLOSEST, + FIND_ELEMENT }; struct gprtRayHit { @@ -49,16 +50,9 @@ class GPRTRayTracer : public RayTracer { // Setup the different shader programs for use with this ray tracer void setup_shaders(); - MeshID find_element(const Position& point) const override - { - fatal_error("Element trees not currently supported with GPRT ray tracer"); - return ID_NONE; - }; + MeshID find_element(const Position& point) const override; - MeshID find_element(TreeID tree, const Position& point) const override { - fatal_error("Element trees not currently supported with GPRT ray tracer"); - return ID_NONE; - }; + MeshID find_element(TreeID tree, const Position& point) const override; std::pair register_volume(const std::shared_ptr& mesh_manager, @@ -114,7 +108,8 @@ class GPRTRayTracer : public RayTracer { // Shader programs std::map> rayGenPrograms_; - GPRTMissOf missProgram_; + GPRTMissOf triangleMissProgram_; + GPRTMissOf tetMissProgram_; GPRTComputeOf aabbTriPopulationProgram_; // aabbTetPopulationProgram_; // tetrahedraGeomType_; //> surface_to_geometry_map_; //ray = gprtBufferGetDevicePointer(rayHitBuffers_.ray); rayGenPIVData->hit = gprtBufferGetDevicePointer(rayHitBuffers_.hit); + dblRayGenData* rayGenFindElementData = gprtRayGenGetParameters(rayGenPrograms_.at(RayGenType::FIND_ELEMENT)); + rayGenFindElementData->ray = gprtBufferGetDevicePointer(rayHitBuffers_.ray); + rayGenFindElementData->hit = gprtBufferGetDevicePointer(rayHitBuffers_.hit); + // Set up build parameters for acceleration structures buildParams_.buildMode = GPRT_BUILD_MODE_FAST_BUILD_NO_UPDATE; } @@ -41,11 +45,16 @@ GPRTRayTracer::~GPRTRayTracer() gprtComputeSynchronize(context_); - // Destroy TLAS structures + // Destroy TLAS structures for surface trees for (const auto& [tree, accel] : surface_volume_tree_to_accel_map) { gprtAccelDestroy(accel); } + // Destroy TLAS structures for element trees + for (const auto& [tree, accel] : element_volume_tree_to_accel_map) { + gprtAccelDestroy(accel); + } + // Destroy BLAS structures for (const auto& blas : blas_handles_) { gprtAccelDestroy(blas); @@ -56,6 +65,7 @@ GPRTRayTracer::~GPRTRayTracer() gprtGeomDestroy(geom); } gprtGeomTypeDestroy(trianglesGeomType_); + gprtGeomTypeDestroy(tetrahedraGeomType_); // Destroy Buffers gprtBufferDestroy(rayHitBuffers_.ray); @@ -72,16 +82,24 @@ void GPRTRayTracer::setup_shaders() // Set up ray generation and miss programs rayGenPrograms_[RayGenType::RAY_FIRE] = gprtRayGenCreate(context_, module_, "ray_fire"); rayGenPrograms_[RayGenType::POINT_IN_VOLUME] = gprtRayGenCreate(context_, module_, "point_in_volume"); + rayGenPrograms_[RayGenType::FIND_ELEMENT] = gprtRayGenCreate(context_, module_, "find_element"); // TODO: Add Occluded and closest raygen entry points - missProgram_ = gprtMissCreate(context_, module_, "ray_fire_miss"); + triangleMissProgram_ = gprtMissCreate(context_, module_, "ray_fire_miss"); + tetMissProgram_ = gprtMissCreate(context_, module_, "tet_miss"); + aabbTriPopulationProgram_ = gprtComputeCreate(context_, module_, "populate_tri_aabbs"); aabbTetPopulationProgram_ = gprtComputeCreate(context_, module_, "populate_tet_aabbs"); // Create a "triangle" geometry type and set its closest-hit program trianglesGeomType_ = gprtGeomTypeCreate(context_, GPRT_AABBS); - gprtGeomTypeSetClosestHitProg(trianglesGeomType_, 0, module_, "ray_fire_hit"); // closesthit for ray queries - gprtGeomTypeSetIntersectionProg(trianglesGeomType_, 0, module_, "DPTrianglePluckerIntersection"); // set intersection program for double precision rays + gprtGeomTypeSetClosestHitProg(trianglesGeomType_, RT_SURFACE_RAY_INDEX, module_, "ray_fire_hit"); // closesthit for ray queries + gprtGeomTypeSetIntersectionProg(trianglesGeomType_, RT_SURFACE_RAY_INDEX, module_, "DPTrianglePluckerIntersection"); // set intersection program for double precision rays against triangles + + // Create a "tetrahedron" geometry type and set its closest-hit program + tetrahedraGeomType_ = gprtGeomTypeCreate(context_, GPRT_AABBS); + gprtGeomTypeSetClosestHitProg(tetrahedraGeomType_, RT_VOLUME_RAY_INDEX, module_, "tet_contain_hit"); // closesthit for point-in-volume queries + gprtGeomTypeSetIntersectionProg(tetrahedraGeomType_, RT_VOLUME_RAY_INDEX, module_, "DPTetrahedronPluckerIntersection"); // set intersection program for double precision rays against tetrahedra } void GPRTRayTracer::init() @@ -173,8 +191,7 @@ GPRTRayTracer::create_surface_tree(const std::shared_ptr& mesh_mana gprtAccelBuild(context_, blas, buildParams_); - gprt::Instance instance; - instance = gprtAccelGetInstance(blas); // create instance of BLAS to be added to TLAS + gprt::Instance instance = gprtAccelGetInstance(blas); // create instance of BLAS to be added to TLAS instance.mask = 0xff; // mask can be used to filter instances during ray traversal. 0xff ensures no filtering // Store in maps @@ -219,7 +236,8 @@ GPRTRayTracer::create_element_tree(const std::shared_ptr& mesh_mana DPTetrahedronGeomData* geom_data = nullptr; auto tetrahedraGeom = gprtGeomCreate(context_, tetrahedraGeomType_); - + geom_data = gprtGeomGetParameters(tetrahedraGeom); // pointer to assign data to + auto vertices = mesh_manager->get_volume_vertices(volume_id); auto indices = mesh_manager->get_volume_connectivity(volume_id); std::vector dbl3Vertices; @@ -235,8 +253,18 @@ GPRTRayTracer::create_element_tree(const std::shared_ptr& mesh_mana ui4Indices.emplace_back(indices[i], indices[i + 1], indices[i + 2], indices[i + 3]); } + // Get storage for prim IDs + std::vector primitive_refs; + primitive_refs.reserve(mesh_manager->num_volume_elements(volume_id)); + for (const auto &element : mesh_manager->get_volume_elements(volume_id)) { + GPRTPrimitiveRef prim_ref; + prim_ref.id = element; + primitive_refs.push_back(prim_ref); + } + auto vertex_buffer = gprtDeviceBufferCreate(context_, dbl3Vertices.size(), dbl3Vertices.data()); auto connectivity_buffer = gprtDeviceBufferCreate(context_, ui4Indices.size(), ui4Indices.data()); + auto primitive_refs_buffer = gprtDeviceBufferCreate(context_, primitive_refs.size(), primitive_refs.data()); auto aabb_buffer = gprtDeviceBufferCreate(context_, 2*volume_elements.size(), 0); // AABBs for each tetrahedron gprtAABBsSetPositions(tetrahedraGeom, aabb_buffer, volume_elements.size(), 2*sizeof(float3), 0); @@ -245,14 +273,21 @@ GPRTRayTracer::create_element_tree(const std::shared_ptr& mesh_mana geom_data->index = gprtBufferGetDevicePointer(connectivity_buffer); geom_data->num_tets = volume_elements.size(); geom_data->vol_id = volume_id; + geom_data->ray = gprtBufferGetDevicePointer(rayHitBuffers_.ray); + geom_data->primitive_refs = gprtBufferGetDevicePointer(primitive_refs_buffer); gprtComputeLaunch(aabbTetPopulationProgram_, {volume_elements.size(), 1, 1}, {1, 1, 1}, *geom_data); - GPRTAccel volume_element_accel = gprtAABBAccelCreate(context_, tetrahedraGeom, buildParams_.buildMode); + GPRTAccel blas = gprtAABBAccelCreate(context_, tetrahedraGeom, buildParams_.buildMode); + gprtAccelBuild(context_, blas, buildParams_); + gprt::Instance instance = gprtAccelGetInstance(blas); // create instance of BLAS to be added to TLAS + instance.mask = 0xff; // mask can be used to filter instances during ray traversal. 0xff ensures no filtering - gprtAccelBuild(context_, volume_element_accel, buildParams_); + auto instanceBuffer = gprtDeviceBufferCreate(context_, 1, &instance); + GPRTAccel volume_tlas = gprtInstanceAccelCreate(context_, 1, instanceBuffer); + gprtAccelBuild(context_, volume_tlas, buildParams_); - // element_volume_tree_to_accel_map[tree] = volume_element_accel; + element_volume_tree_to_accel_map[tree] = volume_tlas; return tree; }; @@ -412,4 +447,39 @@ void GPRTRayTracer::check_ray_buffer_capacity(size_t N) gprtBuildShaderBindingTable(context_, static_cast(GPRT_SBT_GEOM | GPRT_SBT_RAYGEN)); } +MeshID GPRTRayTracer::find_element(const Position& point) const +{ + return find_element(global_element_tree_, point); +} + +MeshID GPRTRayTracer::find_element(TreeID tree, const Position& point) const +{ + + GPRTAccel volume = element_volume_tree_to_accel_map.at(tree); + auto rayGen = rayGenPrograms_.at(RayGenType::FIND_ELEMENT); + dblRayGenData* rayGenData = gprtRayGenGetParameters(rayGen); + + gprtBufferMap(rayHitBuffers_.ray); // Update the ray input buffer + dblRay* ray = gprtBufferGetHostPointer(rayHitBuffers_.ray); + ray[0].volume_accel_solid = gprtAccelGetDeviceAddress(volume); + ray[0].origin = {point.x, point.y, point.z}; + ray[0].volume_tree = tree; // Set the TreeID of the volume being queried + + gprtBufferUnmap(rayHitBuffers_.ray); // required to sync buffer back on GPU? + + gprtRayGenLaunch1D(context_, rayGen, 1); // Launch raygen shader (entry point to RT pipeline) + gprtGraphicsSynchronize(context_); // Ensure all GPU operations are complete before returning control flow to CPU + + // Retrieve the hit from the dblHit buffer + gprtBufferMap(rayHitBuffers_.hit); + dblHit* hit = gprtBufferGetHostPointer(rayHitBuffers_.hit); + auto primitive_id = hit[0].primitive_id; + gprtBufferUnmap(rayHitBuffers_.hit); // required to sync buffer back on GPU? Maybe this second unmap isn't actually needed since we dont need to resyncrhonize after retrieving the data from device + + if (primitive_id == -1) + return ID_NONE; + else + return primitive_id; +} + } // namespace xdg diff --git a/src/gprt/tetrahedron_rt_shaders.slang b/src/gprt/tetrahedron_rt_shaders.slang index d9e9cc95..4224ac72 100644 --- a/src/gprt/tetrahedron_rt_shaders.slang +++ b/src/gprt/tetrahedron_rt_shaders.slang @@ -30,26 +30,26 @@ populate_tet_aabbs(uint3 DispatchThreadID: SV_DispatchThreadID, uniform DPTetrah } [shader("raygeneration")] -void find_element(uniform dblRayGenData record, uniform DPTetrahedronGeomData mesh) { +void find_element(uniform dblRayGenData record) { // This shader is only used to find the element containing a point, so we can ignore the raygen data and just use the ray's origin as the point to test for containment SolidRayFirePayload payload; uint rayID = DispatchRaysIndex().x; PointDesc pointDesc; - pointDesc.Origin = float3(mesh.ray[rayID].origin); + pointDesc.Origin = float3(record.ray[rayID].origin); - SolidAccelerationStructure world = mesh.ray[rayID].volume_accel_solid; + SolidAccelerationStructure world = record.ray[rayID].volume_accel_solid; // Pass the ray's origin and direction to the payload payload.primitiveid = -1; // Initialize to -1 to indicate the point is - TracePoint(world, RAY_FLAG_NONE, 0xff, 0, 1, pointDesc, payload); + TracePoint(world, RAY_FLAG_NONE, 0xff, xdg::RT_VOLUME_RAY_INDEX, xdg::RT_VOLUME_MISS_INDEX, pointDesc, payload); // Store the ID of the containing element in the hit buffer for CPU readback - record.hit[rayID].primitive_id = payload.primitiveid; + record.hit[rayID].primitive_id = payload.primitiveid; // we're reusing dblHit's primitive_id field to store the ID of the containing element since we don't need hit distance or surf_id for this shader } [shader("miss")] -void tet_miss(uniform DPTetrahedronGeomData record, inout SolidRayFirePayload payload) { +void tet_miss(inout SolidRayFirePayload payload) { payload.primitiveid = -1; // Mark the ray as outside all tetrahedra } @@ -79,11 +79,9 @@ void DPTetrahedronPluckerIntersection(uniform DPTetrahedronGeomData record) double3 v1 = record.vertex[indices[1]]; double3 v2 = record.vertex[indices[2]]; double3 v3 = record.vertex[indices[3]]; - double3 origin = record.ray[rayID].origin; bool inside = plucker_tet_containment_test(origin, v0, v1, v2, v3); - if (!inside) return; // If we are inside the tet we need to report a hit to move through the rest of the RT pipeline @@ -92,7 +90,6 @@ void DPTetrahedronPluckerIntersection(uniform DPTetrahedronGeomData record) DPAttribute attr; attr.global_prim_id = global_prim_id; - ReportHit(f32t, hit_kind, attr); } diff --git a/src/gprt/triangle_rt_shaders.slang b/src/gprt/triangle_rt_shaders.slang index 641f4077..6e1727d8 100644 --- a/src/gprt/triangle_rt_shaders.slang +++ b/src/gprt/triangle_rt_shaders.slang @@ -43,7 +43,7 @@ void ray_fire_miss(inout SurfaceRayFirePayload payload) { // This ray generation program will kick off the ray tracing process, // generating rays and tracing them into the world. [shader("raygeneration")] -void ray_fire(uniform dblRayGenData record, uniform DPTriangleGeomData mesh) { +void ray_fire(uniform dblRayGenData record) { SurfaceRayFirePayload payload; uint rayID = DispatchRaysIndex().x; @@ -61,7 +61,7 @@ void ray_fire(uniform dblRayGenData record, uniform DPTriangleGeomData mesh) { payload.surf_id = -1; payload.tlas = world; - TraceRay(world, RAY_FLAG_NONE, 0xff, 0, 1, rayDesc, payload); + TraceRay(world, RAY_FLAG_NONE, 0xff, xdg::RT_SURFACE_RAY_INDEX, xdg::RT_SURFACE_MISS_INDEX, rayDesc, payload); // Store the distance to the hit point and the surface ID in buffers for CPU record.hit[rayID].distance = payload.distance; @@ -70,7 +70,7 @@ void ray_fire(uniform dblRayGenData record, uniform DPTriangleGeomData mesh) { } [shader("raygeneration")] -void point_in_volume(uniform dblRayGenData record, uniform DPTriangleGeomData mesh) { +void point_in_volume(uniform dblRayGenData record) { SurfaceRayFirePayload payload; uint rayID = DispatchRaysIndex().x; @@ -88,7 +88,7 @@ void point_in_volume(uniform dblRayGenData record, uniform DPTriangleGeomData me payload.tlas = world; payload.piv = xdg::PointInVolume::OUTSIDE; // Initialize point in volume check result to outside (0) - TraceRay(world, RAY_FLAG_NONE, 0xff, 0, 1, rayDesc, payload); + TraceRay(world, RAY_FLAG_NONE, 0xff, xdg::RT_SURFACE_RAY_INDEX, xdg::RT_SURFACE_MISS_INDEX, rayDesc, payload); record.hit[rayID].surf_id = payload.surf_id; record.hit[rayID].piv = payload.piv; // Point in volume check result From 63a57c6da9b0e9481e8408c39df25ed2bd77a9c6 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Thu, 19 Feb 2026 14:45:07 +0000 Subject: [PATCH 10/17] Extended find_element test to also test with GPRT backend --- tests/test_find_element.cpp | 73 ++++++++++++++++++++++--------------- 1 file changed, 44 insertions(+), 29 deletions(-) diff --git a/tests/test_find_element.cpp b/tests/test_find_element.cpp index 17023f2b..72ab64b3 100644 --- a/tests/test_find_element.cpp +++ b/tests/test_find_element.cpp @@ -1,45 +1,60 @@ // for testing #include +#include + // xdg includes #include "xdg/constants.h" #include "xdg/mesh_manager_interface.h" #include "xdg/embree/ray_tracer.h" +#include "util.h" #include "mesh_mock.h" using namespace xdg; +using namespace xdg::test; -TEST_CASE("Test Find Volumetric Element") +TEMPLATE_TEST_CASE("Test Find Volumetric Element", "[find_element][mock][volumetric]", + Embree_Raytracer, + GPRT_Raytracer) { - // create a mock mesh manager without volumetric elements - std::shared_ptr mm = std::make_shared(); - mm->init(); // this should do nothing - - REQUIRE(mm->num_volumes() == 1); - REQUIRE(mm->num_surfaces() == 6); - REQUIRE(mm->num_volume_elements(1) == 12); // should return 12 volumetric elements - - std::shared_ptr rti = std::make_shared(); - std::unordered_map volume_to_scene_map; - std::unordered_map element_to_scene_map; - for (auto volume: mm->volumes()) { - auto [volume_tree, element_tree] = rti->register_volume(mm, volume); - volume_to_scene_map[volume] = volume_tree; - element_to_scene_map[volume_tree] = element_tree; + // Generate one test run per enabled backend + constexpr auto rt_backend = TestType::value; + check_ray_tracer_supported(rt_backend); // skip if backend not enabled at configuration time + + DYNAMIC_SECTION(fmt::format("Backend = {}", rt_backend)) + { + auto rti = create_raytracer(rt_backend); + REQUIRE(rti); + // create a mock mesh manager without volumetric elements + std::shared_ptr mm = std::make_shared(); + mm->init(); // this should do nothing + + REQUIRE(mm->num_volumes() == 1); + REQUIRE(mm->num_surfaces() == 6); + REQUIRE(mm->num_volume_elements(1) == 12); // should return 12 volumetric elements + + std::unordered_map volume_to_scene_map; + std::unordered_map element_to_scene_map; + for (auto volume: mm->volumes()) { + auto [volume_tree, element_tree] = rti->register_volume(mm, volume); + volume_to_scene_map[volume] = volume_tree; + element_to_scene_map[volume_tree] = element_tree; + } + REQUIRE(rti->num_registered_trees() == 2); + rti->init(); // Ensure ray tracer is initialized (e.g. build SBT for GPRT after volumes registered) + + MeshID volume = 1; + + // test finding a volumetric element within the volume + Position point_inside {0.0, 0.0, 0.0}; // point inside the volume + MeshID element_id = rti->find_element(element_to_scene_map[volume], point_inside); + REQUIRE(element_id != ID_NONE); // should find an element + REQUIRE(element_id == 7); // Added this hardcode check to ensure GPRT and Embree return the same element ID for the same point + + Position point_outside {10.0, 10.0, 10.0}; // point outside the volume + element_id = rti->find_element(element_to_scene_map[volume], point_outside); + REQUIRE(element_id == ID_NONE); // should not find an element since the point is outside the volume } - REQUIRE(rti->num_registered_trees() == 2); - - MeshID volume = 1; - - // test finding a volumetric element within the volume - Position point_inside {0.0, 0.0, 0.0}; // point inside the volume - MeshID element_id = rti->find_element(element_to_scene_map[volume], point_inside); - REQUIRE(element_id != ID_NONE); // should find an element - - Position point_outside {10.0, 10.0, 10.0}; // point outside the volume - element_id = rti->find_element(element_to_scene_map[volume], point_outside); - REQUIRE(element_id == ID_NONE); // should not find an element since the point is outside the volume } - From 56c29870064074dcd4a28d9b4c5501d3cc9737b1 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Thu, 19 Feb 2026 17:07:39 +0000 Subject: [PATCH 11/17] Added a new slang __generic that populates AABBs for objects of arbitrary numbers of vertices --- include/xdg/gprt/rt_common.slangh | 26 ++++++++++++++++++++++++-- src/gprt/tetrahedron_rt_shaders.slang | 17 ++++------------- src/gprt/triangle_rt_shaders.slang | 16 ++++------------ 3 files changed, 32 insertions(+), 27 deletions(-) diff --git a/include/xdg/gprt/rt_common.slangh b/include/xdg/gprt/rt_common.slangh index 6a37d2ab..4c93f079 100644 --- a/include/xdg/gprt/rt_common.slangh +++ b/include/xdg/gprt/rt_common.slangh @@ -17,10 +17,32 @@ struct DPAttribute static const double PLUCKER_ZERO_TOL = double(20.0) * DBL_EPSILON; -double3x3 inverse3x3(in double3x3 matrix); - // ------------------------------------------------- Helper functions ------------------------------------------------- +// Templated/generic function to compute AABBs in double precision and fill the float AABB buffer. Generic to handle triangles, tets and future primitives with different numbers of vertices +__generic +inline void populate_aabb(uint primID, double3 *vertex, float3 *aabb, uint indices[N]) +{ + double3 dpAabbMin = vertex[indices[0]]; + double3 dpAabbMax = dpAabbMin; + + // unroll the loop since N is known to be small(ish) at compile time to keep performance + [unroll] + for (uint i = 1; i < N; ++i) + { + double3 vert = vertex[indices[i]]; + dpAabbMin = min(dpAabbMin, vert); + dpAabbMax = max(dpAabbMax, vert); + } + + float3 fpAabbMin = float3(dpAabbMin - double3(FLT_EPSILON, FLT_EPSILON, FLT_EPSILON)); + float3 fpAabbMax = float3(dpAabbMax + double3(FLT_EPSILON, FLT_EPSILON, FLT_EPSILON)); + + aabb[2 * primID] = fpAabbMin; + aabb[2 * primID + 1] = fpAabbMax; +} + + bool orientation_cull(in double3 ray, in double3 normal, in xdg::HitOrientation orientation) { if (orientation == xdg::HitOrientation::ANY) return false; // No culling if (orientation == xdg::HitOrientation::EXITING) return dot(ray, normal) < 0.0; // Cull exiting rays diff --git a/src/gprt/tetrahedron_rt_shaders.slang b/src/gprt/tetrahedron_rt_shaders.slang index 4224ac72..0bb177f2 100644 --- a/src/gprt/tetrahedron_rt_shaders.slang +++ b/src/gprt/tetrahedron_rt_shaders.slang @@ -14,19 +14,10 @@ struct SolidRayFirePayload { [numthreads(1, 1, 1)] void populate_tet_aabbs(uint3 DispatchThreadID: SV_DispatchThreadID, uniform DPTetrahedronGeomData record) { - int primID = DispatchThreadID.x; - int4 indices = record.index[primID]; - double3 A = record.vertex[indices[0]]; - double3 B = record.vertex[indices[1]]; - double3 C = record.vertex[indices[2]]; - double3 D = record.vertex[indices[3]]; - double3 dpaabbmin = min(min(A, B), min(C, D)); - double3 dpaabbmax = max(max(A, B), max(C, D)); - float3 fpaabbmin = float3(dpaabbmin - float3(FLT_EPSILON, FLT_EPSILON, FLT_EPSILON)); - float3 fpaabbmax = float3(dpaabbmax + float3(FLT_EPSILON, FLT_EPSILON, FLT_EPSILON)); - - record.aabbs[2 * primID] = fpaabbmin; - record.aabbs[2 * primID + 1] = fpaabbmax; + uint primID = DispatchThreadID.x; + uint4 indices = record.index[primID]; + uint vertexIndices[4] = { indices.x, indices.y, indices.z, indices.w }; + populate_aabb<4>(primID, record.vertex, record.aabbs, vertexIndices); } [shader("raygeneration")] diff --git a/src/gprt/triangle_rt_shaders.slang b/src/gprt/triangle_rt_shaders.slang index 6e1727d8..1c1e434d 100644 --- a/src/gprt/triangle_rt_shaders.slang +++ b/src/gprt/triangle_rt_shaders.slang @@ -100,18 +100,10 @@ void point_in_volume(uniform dblRayGenData record) { [numthreads(1, 1, 1)] void populate_tri_aabbs(uint3 DispatchThreadID: SV_DispatchThreadID, uniform DPTriangleGeomData record) { - int primID = DispatchThreadID.x; - int3 indices = record.index[primID]; - double3 A = record.vertex[indices[0]]; - double3 B = record.vertex[indices[1]]; - double3 C = record.vertex[indices[2]]; - double3 dpaabbmin = min(min(A, B), C); - double3 dpaabbmax = max(max(A, B), C); - float3 fpaabbmin = float3(dpaabbmin - float3(FLT_EPSILON, FLT_EPSILON, FLT_EPSILON)); - float3 fpaabbmax = float3(dpaabbmax + float3(FLT_EPSILON, FLT_EPSILON, FLT_EPSILON)); - - record.aabbs[2 * primID] = fpaabbmin; - record.aabbs[2 * primID + 1] = fpaabbmax; + uint primID = DispatchThreadID.x; + uint3 indices = record.index[primID]; + uint vertexIndices[3] = {indices.x, indices.y, indices.z}; + populate_aabb<3>(primID, record.vertex, record.aabbs, vertexIndices); } // ------------------------------------------------ CUSTOM INTERSECTION SHADERS ------------------------------------------------ From 874c389a92af1fed84085e9b1dafcb22fe2a928c Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Thu, 26 Feb 2026 15:23:18 +0000 Subject: [PATCH 12/17] Changed a constant to be static const to be slang compliant --- include/xdg/geometry/plucker.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/xdg/geometry/plucker.h b/include/xdg/geometry/plucker.h index b35140a4..d08c309d 100644 --- a/include/xdg/geometry/plucker.h +++ b/include/xdg/geometry/plucker.h @@ -25,7 +25,7 @@ struct PluckerIntersectionResult { double t = 0.0; // Distance along the ray to the intersection point }; -constexpr PluckerIntersectionResult EXIT_EARLY = {false, 0.0}; +static const PluckerIntersectionResult EXIT_EARLY = {false, 0.0}; /* Function to return the vertex with the lowest coordinates. To force the same ray-edge computation, the Plücker test needs to use consistent edge From ac2b0f82c3ffb825df8b335c21c56baf10f19eb6 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Fri, 27 Feb 2026 16:38:59 +0000 Subject: [PATCH 13/17] Made plucker_tet_containment code cross-compatible between C++ and slang --- include/xdg/geometry/dp_math.h | 19 ++++---- include/xdg/geometry/plucker.h | 35 +++++++++++++- include/xdg/gprt/rt_common.slangh | 66 --------------------------- include/xdg/tetrahedron_contain.h | 6 +-- src/gprt/tetrahedron_rt_shaders.slang | 2 +- src/tetrahedron_contain.cpp | 36 +-------------- 6 files changed, 46 insertions(+), 118 deletions(-) diff --git a/include/xdg/geometry/dp_math.h b/include/xdg/geometry/dp_math.h index 5b162d3f..1df8b447 100644 --- a/include/xdg/geometry/dp_math.h +++ b/include/xdg/geometry/dp_math.h @@ -8,36 +8,37 @@ functions for dot product, cross product, and absolute value. In C++ compilation while in Slang compilation, it maps to `double3`. */ -#ifdef __SLANG__ +#if defined(__SLANG__) || defined(__SLANG_COMPILER__) // Slang compilation, map dp::vec3 -> double3 namespace dp { typedef double3 vec3; + static const double DBL_EPS = 2.2204460492503131e-016; + static const double PLUCKER_ZERO_TOL = 20.0 * DBL_EPS; + static const double INFTY = 1.7976931348623157e+308; // std::numeric_limits::max() is not available in slang + inline double dot(vec3 a, vec3 b) { return ::dot(a, b); } inline vec3 cross(vec3 a, vec3 b) { return ::cross(a, b); } inline double abs(double a) { return ::abs(a); } - - static const double DBL_ZERO_TOL = 20 * DBL_EPSILON; - static const double INFTY = 1.7976931348623157e+308; // std::numeric_limits::max() is not available in slang } #else #include "xdg/vec3da.h" +#include "xdg/util/linalg.h" // C++ compilation map dp::vec3 -> xdg::Vec3da namespace dp { typedef xdg::Vec3da vec3; + static constexpr double PLUCKER_ZERO_TOL = 20.0 * std::numeric_limits::epsilon(); + constexpr double INFTY {std::numeric_limits::max()}; + inline double dot(vec3 a, vec3 b) { return xdg::dot(a, b); } inline vec3 cross(vec3 a, vec3 b) { return xdg::cross(a, b); } inline double abs(double a) { return std::fabs(a); } - - static constexpr double DBL_ZERO_TOL = 20.0 * std::numeric_limits::epsilon(); - constexpr double INFTY {std::numeric_limits::max()}; - } #endif // end of ifdef __slang__ -#endif // DP_MATH_H \ No newline at end of file +#endif // DP_MATH_H diff --git a/include/xdg/geometry/plucker.h b/include/xdg/geometry/plucker.h index d08c309d..200bb4db 100644 --- a/include/xdg/geometry/plucker.h +++ b/include/xdg/geometry/plucker.h @@ -56,7 +56,7 @@ inline double plucker_edge_test(dp::vec3 vertexa, dp::vec3 vertexb, pip = dp::dot(ray, edge_normal) + dp::dot(ray_normal, edge); pip = -pip; } - if (dp::DBL_ZERO_TOL > dp::abs(pip)) // <-- absd + if (dp::PLUCKER_ZERO_TOL > dp::abs(pip)) // <-- absd pip = 0.0; return pip; } @@ -160,6 +160,37 @@ inline PluckerIntersectionResult plucker_ray_tri_intersect(dp::vec3 vertices[3], return {true, dist_out}; } +// Plücker containment test for tetrahedra. Returns true if the point is inside the tetrahedron defined by vertices v0, v1, v2, v3. +inline bool plucker_tet_containment_test(const dp::vec3 point, + const dp::vec3 v0, + const dp::vec3 v1, + const dp::vec3 v2, + const dp::vec3 v3) { + + // TODO - I've decided to use Cramer's rule here instead of matrix inversion to make it easier to implement a cross-compilable version of this function + // Not sure if that is necessarily the best choice however. + const dp::vec3 e0 = v1 - v0; + const dp::vec3 e1 = v2 - v0; + const dp::vec3 e2 = v3 - v0; + const dp::vec3 rhs = point - v0; + const double det = dp::dot(e0, dp::cross(e1, e2)); // scalar triple product of matrix [e0 e1 e2] + + const double inv_det = 1.0 / det; + const double lambda1 = dp::dot(rhs, dp::cross(e1, e2)) * inv_det; + const double lambda2 = dp::dot(e0, dp::cross(rhs, e2)) * inv_det; + const double lambda3 = dp::dot(e0, dp::cross(e1, rhs)) * inv_det; + const double lambda0 = 1.0 - (lambda1 + lambda2 + lambda3); + + const double barycentric_min = -dp::PLUCKER_ZERO_TOL; + const double barycentric_max = 1.0 + dp::PLUCKER_ZERO_TOL; + + // Check all λ_i in [0, 1] + return (lambda0 >= barycentric_min && lambda0 <= barycentric_max) && + (lambda1 >= barycentric_min && lambda1 <= barycentric_max) && + (lambda2 >= barycentric_min && lambda2 <= barycentric_max) && + (lambda3 >= barycentric_min && lambda3 <= barycentric_max); +} + } // namespace xdg -#endif // include guard \ No newline at end of file +#endif // include guard diff --git a/include/xdg/gprt/rt_common.slangh b/include/xdg/gprt/rt_common.slangh index 4c93f079..7f93ded3 100644 --- a/include/xdg/gprt/rt_common.slangh +++ b/include/xdg/gprt/rt_common.slangh @@ -15,8 +15,6 @@ struct DPAttribute int global_prim_id; }; -static const double PLUCKER_ZERO_TOL = double(20.0) * DBL_EPSILON; - // ------------------------------------------------- Helper functions ------------------------------------------------- // Templated/generic function to compute AABBs in double precision and fill the float AABB buffer. Generic to handle triangles, tets and future primitives with different numbers of vertices @@ -76,68 +74,4 @@ float next_after(float a) { return asfloat(a_); } -// slang doesnt seem to have an inverse matrix intrinsic so I needed to write my own unfortunately :/ -double3x3 inverse3x3(in double3x3 matrix) -{ - double m00 = matrix[0][0], m01 = matrix[0][1], m02 = matrix[0][2]; - double m10 = matrix[1][0], m11 = matrix[1][1], m12 = matrix[1][2]; - double m20 = matrix[2][0], m21 = matrix[2][1], m22 = matrix[2][2]; - - // Cofactor matrix. - double c00 = (m11 * m22 - m12 * m21); - double c01 = -(m10 * m22 - m12 * m20); - double c02 = (m10 * m21 - m11 * m20); - double c10 = -(m01 * m22 - m02 * m21); - double c11 = (m00 * m22 - m02 * m20); - double c12 = -(m00 * m21 - m01 * m20); - double c20 = (m01 * m12 - m02 * m11); - double c21 = -(m00 * m12 - m02 * m10); - double c22 = (m00 * m11 - m01 * m10); - - double det = m00 * c00 + m01 * c01 + m02 * c02; - - // Adjugate is transpose(cofactor). - double3x3 adj = { - { c00, c10, c20 }, - { c01, c11, c21 }, - { c02, c12, c22 } - }; - return adj * (double(1.0) / det); -} - -bool plucker_tet_containment_test(in double3 point, - in double3 v0, - in double3 v1, - in double3 v2, - in double3 v3) { - // Create matrix T = [v1 - v0, v2 - v0, v3 - v0] - double3 e0 = v1 - v0; - double3 e1 = v2 - v0; - double3 e2 = v3 - v0; - double3x3 T = { { e0.x, e0.y, e0.z }, - { e1.x, e1.y, e1.z }, - { e2.x, e2.y, e2.z } }; - - // Vector from v0 to point - double3 rhs = point - v0; - - // Solve T * [λ1, λ2, λ3] = rhs. - // In Slang, mul(matrix, vector) is lowered as vector-times-matrix on shader backends, - // so transpose the inverse to preserve CPU-side convention. - double3 lambda123 = mul(transpose(inverse3x3(T)), double3(rhs.x, rhs.y, rhs.z)); - - // Compute λ0 - double lambda0 = double(1.0) - (lambda123.x + lambda123.y + lambda123.z); - - // Final barycentric coordinate vector - double4 bary = { lambda0, lambda123.x, lambda123.y, lambda123.z }; - - // Check all λ_i in [0, 1] - for (int i = 0; i < 4; ++i) { - if (bary[i] < -PLUCKER_ZERO_TOL || bary[i] > double(1.0) + PLUCKER_ZERO_TOL) - return false; - } - return true; -} - #endif // XDG_GPRT_RT_COMMON_SLANGH diff --git a/include/xdg/tetrahedron_contain.h b/include/xdg/tetrahedron_contain.h index 51f1a338..03250447 100644 --- a/include/xdg/tetrahedron_contain.h +++ b/include/xdg/tetrahedron_contain.h @@ -3,6 +3,7 @@ #include "xdg/vec3da.h" +#include "xdg/geometry/plucker.h" namespace xdg { @@ -30,11 +31,6 @@ namespace xdg * @return `true` if the point is inside or on the boundary of the tetrahedron, * `false` otherwise. */ -bool plucker_tet_containment_test(const Position& point, - const Vertex& v0, - const Vertex& v1, - const Vertex& v2, - const Vertex& v3); // Embree call back functions for element search void VolumeElementBoundsFunc(RTCBoundsFunctionArguments* args); diff --git a/src/gprt/tetrahedron_rt_shaders.slang b/src/gprt/tetrahedron_rt_shaders.slang index 0bb177f2..1ecced56 100644 --- a/src/gprt/tetrahedron_rt_shaders.slang +++ b/src/gprt/tetrahedron_rt_shaders.slang @@ -72,7 +72,7 @@ void DPTetrahedronPluckerIntersection(uniform DPTetrahedronGeomData record) double3 v3 = record.vertex[indices[3]]; double3 origin = record.ray[rayID].origin; - bool inside = plucker_tet_containment_test(origin, v0, v1, v2, v3); + bool inside = xdg::plucker_tet_containment_test(origin, v0, v1, v2, v3); if (!inside) return; // If we are inside the tet we need to report a hit to move through the rest of the RT pipeline diff --git a/src/tetrahedron_contain.cpp b/src/tetrahedron_contain.cpp index 743ea1d0..95dfc0b3 100644 --- a/src/tetrahedron_contain.cpp +++ b/src/tetrahedron_contain.cpp @@ -2,46 +2,12 @@ #include "xdg/ray_tracing_interface.h" #include "xdg/ray.h" #include "xdg/vec3da.h" - +#include "xdg/geometry/plucker.h" #include "xdg/util/linalg.h" namespace xdg { -bool plucker_tet_containment_test(const Position& point, - const Position& v0, - const Position& v1, - const Position& v2, - const Position& v3) { - using namespace linalg::aliases; - // Create matrix T = [v1 - v0, v2 - v0, v3 - v0] - Vec3da e0 = v1 - v0; - Vec3da e1 = v2 - v0; - Vec3da e2 = v3 - v0; - double3x3 T = { {e0.x, e0.y, e0.z}, - {e1.x, e1.y, e1.z}, - {e2.x, e2.y, e2.z}}; - - // Vector from v0 to point - Vec3da rhs = point - v0; - - // Solve T * [λ1, λ2, λ3] = rhs - double3 lambda123 = mul(inverse(T),{rhs.x, rhs.y, rhs.z}); - - // Compute λ0 - double lambda0 = 1.0f - (lambda123.x + lambda123.y + lambda123.z); - - // Final barycentric coordinate vector - double4 bary = { lambda0, lambda123.x, lambda123.y, lambda123.z }; - - // Check all λ_i in [0, 1] - for (int i = 0; i < 4; ++i) { - if (bary[i] < -PLUCKER_ZERO_TOL || bary[i] > 1.0f + PLUCKER_ZERO_TOL) - return false; - } - return true; -} - // Embree callbacks void VolumeElementBoundsFunc(RTCBoundsFunctionArguments* args) From 342d07e245e2ca42a52ccabc05a196fa9c9d42f5 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Fri, 27 Feb 2026 17:47:15 +0000 Subject: [PATCH 14/17] Switch to pre-allocate arrays in BVH construction rather than relying on push_back --- include/xdg/geometry/dp_math.h | 1 - src/gprt/ray_tracer.cpp | 72 ++++++++++++++-------------------- 2 files changed, 30 insertions(+), 43 deletions(-) diff --git a/include/xdg/geometry/dp_math.h b/include/xdg/geometry/dp_math.h index 1df8b447..9820934c 100644 --- a/include/xdg/geometry/dp_math.h +++ b/include/xdg/geometry/dp_math.h @@ -25,7 +25,6 @@ namespace dp { #else #include "xdg/vec3da.h" -#include "xdg/util/linalg.h" // C++ compilation map dp::vec3 -> xdg::Vec3da namespace dp { diff --git a/src/gprt/ray_tracer.cpp b/src/gprt/ray_tracer.cpp index ce5e7835..e3860528 100644 --- a/src/gprt/ray_tracer.cpp +++ b/src/gprt/ray_tracer.cpp @@ -128,14 +128,6 @@ GPRTRayTracer::create_surface_tree(const std::shared_ptr& mesh_mana std::vector surfaceBlasInstances; // BLAS for each (surface) geometry in this volume for (const auto &surf : volume_surfaces) { - auto num_faces = mesh_manager->num_surface_faces(surf); - - // get the sense of this surface with respect to the volume - Sense triangle_sense {Sense::UNSET}; - auto surf_to_vol_senses = mesh_manager->get_parent_volumes(surf); - if (volume_id == surf_to_vol_senses.first) triangle_sense = Sense::FORWARD; - else if (volume_id == surf_to_vol_senses.second) triangle_sense = Sense::REVERSE; - DPTriangleGeomData* geom_data = nullptr; auto triangleGeom = gprtGeomCreate(context_, trianglesGeomType_); geom_data = gprtGeomGetParameters(triangleGeom); // pointer to assign data to @@ -143,30 +135,32 @@ GPRTRayTracer::create_surface_tree(const std::shared_ptr& mesh_mana // Get storage for vertices and indices auto vertices = mesh_manager->get_surface_vertices(surf); auto indices = mesh_manager->get_surface_connectivity(surf); - std::vector dbl3Vertices; - dbl3Vertices.reserve(vertices.size()); - for (const auto &vertex : vertices) { - dbl3Vertices.push_back({vertex.x, vertex.y, vertex.z}); + if (indices.size() % 3 != 0) { + fatal_error("Surface {} connectivity size ({}) is not divisible by 3 for triangles", surf, indices.size()); + } + std::vector dbl3Vertices(vertices.size()); + for (size_t i = 0; i < vertices.size(); ++i) { + const auto& vertex = vertices[i]; + dbl3Vertices[i] = {vertex.x, vertex.y, vertex.z}; } // Get storage for indices - std::vector ui3Indices; - ui3Indices.reserve(indices.size() / 3); - for (size_t i = 0; i < indices.size(); i += 3) { - ui3Indices.emplace_back(indices[i], indices[i + 1], indices[i + 2]); + std::vector ui3Indices(indices.size() / 3); + for (size_t i = 0, triIdx = 0; i < indices.size(); i += 3, ++triIdx) { + ui3Indices[triIdx] = uint3(indices[i], indices[i + 1], indices[i + 2]); } // Get storage for normals - std::vector normals; - std::vector primitive_refs; - primitive_refs.reserve(num_faces); - normals.reserve(num_faces); - for (const auto &face : mesh_manager->get_surface_faces(surf)) { + auto surface_faces = mesh_manager->get_surface_faces(surf); + const size_t num_faces = surface_faces.size(); + + std::vector normals(surface_faces.size()); + std::vector primitive_refs(surface_faces.size()); + for (size_t i = 0; i < surface_faces.size(); ++i) { + const auto face = surface_faces[i]; auto norm = mesh_manager->face_normal(face); - normals.push_back({norm.x, norm.y, norm.z}); - GPRTPrimitiveRef prim_ref; - prim_ref.id = face; - primitive_refs.push_back(prim_ref); + normals[i] = {norm.x, norm.y, norm.z}; + primitive_refs[i].id = face; } auto vertex_buffer = gprtDeviceBufferCreate(context_, dbl3Vertices.size(), dbl3Vertices.data()); @@ -197,9 +191,6 @@ GPRTRayTracer::create_surface_tree(const std::shared_ptr& mesh_mana // Store in maps surface_to_geometry_map_[surf] = triangleGeom; - geom_data = gprtGeomGetParameters(triangleGeom); - instance = gprtAccelGetInstance(blas); - instance.mask = 0xff; surfaceBlasInstances.push_back(instance); globalBlasInstances_.push_back(instance); @@ -240,26 +231,23 @@ GPRTRayTracer::create_element_tree(const std::shared_ptr& mesh_mana auto vertices = mesh_manager->get_volume_vertices(volume_id); auto indices = mesh_manager->get_volume_connectivity(volume_id); - std::vector dbl3Vertices; - dbl3Vertices.reserve(vertices.size()); - for (const auto &vertex : vertices) { - dbl3Vertices.push_back({vertex.x, vertex.y, vertex.z}); + + std::vector dbl3Vertices(vertices.size()); + for (size_t i = 0; i < vertices.size(); ++i) { + const auto& vertex = vertices[i]; + dbl3Vertices[i] = {vertex.x, vertex.y, vertex.z}; } // Get storage for indices - std::vector ui4Indices; - ui4Indices.reserve(indices.size() / 4); - for (size_t i = 0; i + 3 < indices.size(); i += 4) { - ui4Indices.emplace_back(indices[i], indices[i + 1], indices[i + 2], indices[i + 3]); + std::vector ui4Indices(indices.size() / 4); + for (size_t i = 0, tetIdx = 0; i < indices.size(); i += 4, ++tetIdx) { + ui4Indices[tetIdx] = uint4(indices[i], indices[i + 1], indices[i + 2], indices[i + 3]); } // Get storage for prim IDs - std::vector primitive_refs; - primitive_refs.reserve(mesh_manager->num_volume_elements(volume_id)); - for (const auto &element : mesh_manager->get_volume_elements(volume_id)) { - GPRTPrimitiveRef prim_ref; - prim_ref.id = element; - primitive_refs.push_back(prim_ref); + std::vector primitive_refs(volume_elements.size()); + for (size_t i = 0; i < volume_elements.size(); ++i) { + primitive_refs[i].id = volume_elements[i]; } auto vertex_buffer = gprtDeviceBufferCreate(context_, dbl3Vertices.size(), dbl3Vertices.data()); From ac4c7f32ee3fd05d0a13f44d386f888fb1fafe17 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Thu, 19 Mar 2026 09:17:45 +0000 Subject: [PATCH 15/17] Rebased to latest main --- src/mesh_manager_interface.cpp | 2 +- vendor/GPRT | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/mesh_manager_interface.cpp b/src/mesh_manager_interface.cpp index 1f192569..82cba68b 100644 --- a/src/mesh_manager_interface.cpp +++ b/src/mesh_manager_interface.cpp @@ -299,7 +299,7 @@ std::vector MeshManager::get_surface_vertices(MeshID surface) const return surface_local_mesh_data(surface).vertices; } -std::vector MeshManager::get_volume_connectivity(MeshID volume) const +std::vector MeshManager::get_surface_connectivity(MeshID surface) const { return surface_local_mesh_data(surface).connectivity; } diff --git a/vendor/GPRT b/vendor/GPRT index f1e95e41..405d9ee9 160000 --- a/vendor/GPRT +++ b/vendor/GPRT @@ -1 +1 @@ -Subproject commit f1e95e4188cde591547d6b4a33a70bf2afaeec59 +Subproject commit 405d9ee9f5ee8e1a0455f776f9e2c3adffb64160 From bb61bf667a625b0e50eebd804284da9dfd6e0449 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Fri, 10 Apr 2026 11:31:25 +0100 Subject: [PATCH 16/17] Addressed comments and created new shared_constants.h header --- CMakeLists.txt | 2 +- include/xdg/constants.h | 18 +------------- include/xdg/geometry/dp_math.h | 2 +- include/xdg/gprt/ray_tracer.h | 2 +- include/xdg/gprt/rt_common.slangh | 5 ++-- include/xdg/gprt/shared_structs.h | 2 +- include/xdg/shared_constants.h | 35 +++++++++++++++++++++++++++ include/xdg/shared_enums.h | 20 --------------- src/gprt/ray_tracer.cpp | 5 +--- src/gprt/tetrahedron_rt_shaders.slang | 10 ++++---- src/gprt/triangle_rt_shaders.slang | 13 ++++------ 11 files changed, 54 insertions(+), 60 deletions(-) create mode 100644 include/xdg/shared_constants.h delete mode 100644 include/xdg/shared_enums.h diff --git a/CMakeLists.txt b/CMakeLists.txt index a9f1afc1..b61e423a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -154,7 +154,7 @@ src/gprt/triangle_rt_shaders.slang src/gprt/tetrahedron_rt_shaders.slang include/xdg/gprt/shared_structs.h include/xdg/gprt/rt_common.slangh -include/xdg/shared_enums.h +include/xdg/shared_constants.h include/xdg/geometry/dp_math.h ) endif() diff --git a/include/xdg/constants.h b/include/xdg/constants.h index 38df9db1..a0052737 100644 --- a/include/xdg/constants.h +++ b/include/xdg/constants.h @@ -8,7 +8,7 @@ #include "fmt/format.h" -#include "xdg/shared_enums.h" +#include "xdg/shared_constants.h" namespace xdg { @@ -70,22 +70,6 @@ static const std::map RT_LIB_TO_STR = {RTLibrary::GPRT, "GPRT"} }; -// Mesh identifer type -using MeshID = int32_t; -using MeshIndex = int32_t; - -// Null mesh ID -constexpr MeshID ID_NONE {-1}; -constexpr MeshIndex INDEX_NONE {-1}; - -// Scene/Tree ID -using TreeID = int32_t; -using SurfaceTreeID = TreeID; -using ElementTreeID = TreeID; - -// Null tree ID -constexpr TreeID TREE_NONE {-1}; - // for abs(x) >= min_rcp_input the newton raphson rcp calculation does not fail constexpr float min_rcp_input = std::numeric_limits::min() /* FIX ME */ *1E5 /* SHOULDNT NEED TO MULTIPLY BY THIS VALUE */; constexpr int BVH_MAX_DEPTH = 64; diff --git a/include/xdg/geometry/dp_math.h b/include/xdg/geometry/dp_math.h index 9820934c..a2fece8f 100644 --- a/include/xdg/geometry/dp_math.h +++ b/include/xdg/geometry/dp_math.h @@ -16,7 +16,7 @@ namespace dp { static const double DBL_EPS = 2.2204460492503131e-016; static const double PLUCKER_ZERO_TOL = 20.0 * DBL_EPS; - static const double INFTY = 1.7976931348623157e+308; // std::numeric_limits::max() is not available in slang + static const double INFTY = 1.7976931348623157e+308; inline double dot(vec3 a, vec3 b) { return ::dot(a, b); } inline vec3 cross(vec3 a, vec3 b) { return ::cross(a, b); } diff --git a/include/xdg/gprt/ray_tracer.h b/include/xdg/gprt/ray_tracer.h index 4e759fe6..13ad8956 100644 --- a/include/xdg/gprt/ray_tracer.h +++ b/include/xdg/gprt/ray_tracer.h @@ -123,7 +123,7 @@ class GPRTRayTracer : public RayTracer { GPRTGeomTypeOf tetrahedraGeomType_; //> surface_to_geometry_map_; // 0.0; // Cull entering rays + double dot_product = dot(ray, normal); + if (orientation == xdg::HitOrientation::EXITING) return dot_product < 0.0; // Cull exiting rays + if (orientation == xdg::HitOrientation::ENTERING) return dot_product > 0.0; // Cull entering rays return false; // Default case, no culling } diff --git a/include/xdg/gprt/shared_structs.h b/include/xdg/gprt/shared_structs.h index 3ec4cba4..7871b939 100644 --- a/include/xdg/gprt/shared_structs.h +++ b/include/xdg/gprt/shared_structs.h @@ -1,5 +1,5 @@ #include "gprt.h" -#include "../shared_enums.h" +#include "../shared_constants.h" #include "../geometry/dp_math.h" diff --git a/include/xdg/shared_constants.h b/include/xdg/shared_constants.h new file mode 100644 index 00000000..53d28706 --- /dev/null +++ b/include/xdg/shared_constants.h @@ -0,0 +1,35 @@ +#ifndef XDG_SHARED_CONSTANTS_H +#define XDG_SHARED_CONSTANTS_H + +namespace xdg { + + // Mesh identifer type + typedef int32_t MeshID; + typedef int32_t MeshIndex; + + // Null mesh ID + static const MeshID ID_NONE = -1; + static const MeshIndex INDEX_NONE = -1; + + // Scene/Tree ID + typedef int32_t TreeID; + typedef TreeID SurfaceTreeID; + typedef TreeID ElementTreeID; + + static const TreeID TREE_NONE = -1; + + enum PointInVolume : int { + UNSET = -1, + OUTSIDE = 0, + INSIDE = 1 + }; + + enum HitOrientation : int { + ANY = -1, + EXITING = 0, + ENTERING = 1, + }; + +} + +#endif // XDG_SHARED_CONSTANTS_H diff --git a/include/xdg/shared_enums.h b/include/xdg/shared_enums.h deleted file mode 100644 index a20c27fa..00000000 --- a/include/xdg/shared_enums.h +++ /dev/null @@ -1,20 +0,0 @@ -#ifndef XDG_SHARED_ENUMS_H -#define XDG_SHARED_ENUMS_H - -namespace xdg { - - enum PointInVolume : int { - UNSET = -1, - OUTSIDE = 0, - INSIDE = 1 - }; - - enum HitOrientation : int { - ANY = -1, - EXITING = 0, - ENTERING = 1, - }; - -} - -#endif // XDG_SHARED_ENUMS_H \ No newline at end of file diff --git a/src/gprt/ray_tracer.cpp b/src/gprt/ray_tracer.cpp index e3860528..6b6db475 100644 --- a/src/gprt/ray_tracer.cpp +++ b/src/gprt/ray_tracer.cpp @@ -464,10 +464,7 @@ MeshID GPRTRayTracer::find_element(TreeID tree, const Position& point) const auto primitive_id = hit[0].primitive_id; gprtBufferUnmap(rayHitBuffers_.hit); // required to sync buffer back on GPU? Maybe this second unmap isn't actually needed since we dont need to resyncrhonize after retrieving the data from device - if (primitive_id == -1) - return ID_NONE; - else - return primitive_id; + return primitive_id; } } // namespace xdg diff --git a/src/gprt/tetrahedron_rt_shaders.slang b/src/gprt/tetrahedron_rt_shaders.slang index 1ecced56..05e512a8 100644 --- a/src/gprt/tetrahedron_rt_shaders.slang +++ b/src/gprt/tetrahedron_rt_shaders.slang @@ -4,7 +4,7 @@ #include "../../include/xdg/gprt/rt_common.slangh" struct SolidRayFirePayload { - int primitiveid; // ID of the primitive the point is contained within + int primitive_id; // ID of the primitive the point is contained within }; @@ -31,22 +31,22 @@ void find_element(uniform dblRayGenData record) { SolidAccelerationStructure world = record.ray[rayID].volume_accel_solid; // Pass the ray's origin and direction to the payload - payload.primitiveid = -1; // Initialize to -1 to indicate the point is + payload.primitive_id = xdg::ID_NONE; // Initialize to ID_NONE to indicate the point is outside all elements TracePoint(world, RAY_FLAG_NONE, 0xff, xdg::RT_VOLUME_RAY_INDEX, xdg::RT_VOLUME_MISS_INDEX, pointDesc, payload); // Store the ID of the containing element in the hit buffer for CPU readback - record.hit[rayID].primitive_id = payload.primitiveid; // we're reusing dblHit's primitive_id field to store the ID of the containing element since we don't need hit distance or surf_id for this shader + record.hit[rayID].primitive_id = payload.primitive_id; // we're reusing dblHit's primitive_id field to store the ID of the containing element since we don't need hit distance or surf_id for this shader } [shader("miss")] void tet_miss(inout SolidRayFirePayload payload) { - payload.primitiveid = -1; // Mark the ray as outside all tetrahedra + payload.primitive_id = xdg::ID_NONE; // Mark the ray as outside all tetrahedra } [shader("closesthit")] void tet_contain_hit(uniform DPTetrahedronGeomData record, inout SolidRayFirePayload payload, in DPAttribute attr) { - payload.primitiveid = attr.global_prim_id; // Store the ID of the tet that contains the point + payload.primitive_id = attr.global_prim_id; // Store the ID of the tet that contains the point } diff --git a/src/gprt/triangle_rt_shaders.slang b/src/gprt/triangle_rt_shaders.slang index 1c1e434d..600f6411 100644 --- a/src/gprt/triangle_rt_shaders.slang +++ b/src/gprt/triangle_rt_shaders.slang @@ -17,9 +17,6 @@ void ray_fire_hit(uniform DPTriangleGeomData record, inout SurfaceRayFirePayload uint hit_kind = HitKind(); uint rayID = DispatchRaysIndex().x; - // There is some logic for handling next volumes inside the h5m-reader which I could make use of too - // TODO : Should the dblHit struct return the next volume ID for the ray back to the host - payload.piv = (hit_kind == HIT_KIND_TRIANGLE_FRONT_FACE) ? xdg::PointInVolume::OUTSIDE : xdg::PointInVolume::INSIDE; @@ -35,8 +32,8 @@ void ray_fire_hit(uniform DPTriangleGeomData record, inout SurfaceRayFirePayload void ray_fire_miss(inout SurfaceRayFirePayload payload) { // Set the miss payload to default values payload.distance = -1.0f; - payload.surf_id = -1; - payload.primitive_id = -1; + payload.surf_id = xdg::ID_NONE; + payload.primitive_id = xdg::ID_NONE; payload.piv = xdg::PointInVolume::UNSET; } @@ -58,7 +55,7 @@ void ray_fire(uniform dblRayGenData record) { // Pass the ray's origin and direction to the payload payload.distance = -1.0f; - payload.surf_id = -1; + payload.surf_id = xdg::ID_NONE; payload.tlas = world; TraceRay(world, RAY_FLAG_NONE, 0xff, xdg::RT_SURFACE_RAY_INDEX, xdg::RT_SURFACE_MISS_INDEX, rayDesc, payload); @@ -84,7 +81,7 @@ void point_in_volume(uniform dblRayGenData record) { SurfaceAccelerationStructure world = record.ray[rayID].volume_accel_surf; // Pass the ray's origin and direction to the payload - payload.surf_id = -1; + payload.surf_id = xdg::ID_NONE; payload.tlas = world; payload.piv = xdg::PointInVolume::OUTSIDE; // Initialize point in volume check result to outside (0) @@ -108,7 +105,7 @@ populate_tri_aabbs(uint3 DispatchThreadID: SV_DispatchThreadID, uniform DPTriang // ------------------------------------------------ CUSTOM INTERSECTION SHADERS ------------------------------------------------ -/* 1D ray generation intersection with a double precision triangle using the Plucker intersection algorithm*/ +// Custom FP64 Plucker ray-triangle intersection algorithm for each ray-triangle pair [shader("intersection")] void DPTrianglePluckerIntersection(uniform DPTriangleGeomData record) { From ace02e9743bcc82f1eef5685fb4211fce9fb70d3 Mon Sep 17 00:00:00 2001 From: waqar-ukaea Date: Mon, 13 Apr 2026 11:11:31 +0100 Subject: [PATCH 17/17] Reverted change to constants --- CMakeLists.txt | 2 +- include/xdg/constants.h | 18 +++++++++++++- include/xdg/gprt/rt_common.slangh | 1 + include/xdg/gprt/shared_structs.h | 2 +- include/xdg/shared_constants.h | 35 --------------------------- include/xdg/shared_enums.h | 20 +++++++++++++++ src/gprt/tetrahedron_rt_shaders.slang | 4 +-- src/gprt/triangle_rt_shaders.slang | 8 +++--- vendor/GPRT | 2 +- 9 files changed, 47 insertions(+), 45 deletions(-) delete mode 100644 include/xdg/shared_constants.h create mode 100644 include/xdg/shared_enums.h diff --git a/CMakeLists.txt b/CMakeLists.txt index b61e423a..a9f1afc1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -154,7 +154,7 @@ src/gprt/triangle_rt_shaders.slang src/gprt/tetrahedron_rt_shaders.slang include/xdg/gprt/shared_structs.h include/xdg/gprt/rt_common.slangh -include/xdg/shared_constants.h +include/xdg/shared_enums.h include/xdg/geometry/dp_math.h ) endif() diff --git a/include/xdg/constants.h b/include/xdg/constants.h index a0052737..38df9db1 100644 --- a/include/xdg/constants.h +++ b/include/xdg/constants.h @@ -8,7 +8,7 @@ #include "fmt/format.h" -#include "xdg/shared_constants.h" +#include "xdg/shared_enums.h" namespace xdg { @@ -70,6 +70,22 @@ static const std::map RT_LIB_TO_STR = {RTLibrary::GPRT, "GPRT"} }; +// Mesh identifer type +using MeshID = int32_t; +using MeshIndex = int32_t; + +// Null mesh ID +constexpr MeshID ID_NONE {-1}; +constexpr MeshIndex INDEX_NONE {-1}; + +// Scene/Tree ID +using TreeID = int32_t; +using SurfaceTreeID = TreeID; +using ElementTreeID = TreeID; + +// Null tree ID +constexpr TreeID TREE_NONE {-1}; + // for abs(x) >= min_rcp_input the newton raphson rcp calculation does not fail constexpr float min_rcp_input = std::numeric_limits::min() /* FIX ME */ *1E5 /* SHOULDNT NEED TO MULTIPLY BY THIS VALUE */; constexpr int BVH_MAX_DEPTH = 64; diff --git a/include/xdg/gprt/rt_common.slangh b/include/xdg/gprt/rt_common.slangh index 5210daa3..cddec57d 100644 --- a/include/xdg/gprt/rt_common.slangh +++ b/include/xdg/gprt/rt_common.slangh @@ -9,6 +9,7 @@ [[vk::push_constant]] dblRayFirePushConstants PC; +static const int32_t ID_NONE = -1; // slang copy of the C++ constant for invalid IDs struct DPAttribute { double f64t; // double precision hit distance diff --git a/include/xdg/gprt/shared_structs.h b/include/xdg/gprt/shared_structs.h index 7871b939..3ec4cba4 100644 --- a/include/xdg/gprt/shared_structs.h +++ b/include/xdg/gprt/shared_structs.h @@ -1,5 +1,5 @@ #include "gprt.h" -#include "../shared_constants.h" +#include "../shared_enums.h" #include "../geometry/dp_math.h" diff --git a/include/xdg/shared_constants.h b/include/xdg/shared_constants.h deleted file mode 100644 index 53d28706..00000000 --- a/include/xdg/shared_constants.h +++ /dev/null @@ -1,35 +0,0 @@ -#ifndef XDG_SHARED_CONSTANTS_H -#define XDG_SHARED_CONSTANTS_H - -namespace xdg { - - // Mesh identifer type - typedef int32_t MeshID; - typedef int32_t MeshIndex; - - // Null mesh ID - static const MeshID ID_NONE = -1; - static const MeshIndex INDEX_NONE = -1; - - // Scene/Tree ID - typedef int32_t TreeID; - typedef TreeID SurfaceTreeID; - typedef TreeID ElementTreeID; - - static const TreeID TREE_NONE = -1; - - enum PointInVolume : int { - UNSET = -1, - OUTSIDE = 0, - INSIDE = 1 - }; - - enum HitOrientation : int { - ANY = -1, - EXITING = 0, - ENTERING = 1, - }; - -} - -#endif // XDG_SHARED_CONSTANTS_H diff --git a/include/xdg/shared_enums.h b/include/xdg/shared_enums.h new file mode 100644 index 00000000..a20c27fa --- /dev/null +++ b/include/xdg/shared_enums.h @@ -0,0 +1,20 @@ +#ifndef XDG_SHARED_ENUMS_H +#define XDG_SHARED_ENUMS_H + +namespace xdg { + + enum PointInVolume : int { + UNSET = -1, + OUTSIDE = 0, + INSIDE = 1 + }; + + enum HitOrientation : int { + ANY = -1, + EXITING = 0, + ENTERING = 1, + }; + +} + +#endif // XDG_SHARED_ENUMS_H \ No newline at end of file diff --git a/src/gprt/tetrahedron_rt_shaders.slang b/src/gprt/tetrahedron_rt_shaders.slang index 05e512a8..47df4479 100644 --- a/src/gprt/tetrahedron_rt_shaders.slang +++ b/src/gprt/tetrahedron_rt_shaders.slang @@ -31,7 +31,7 @@ void find_element(uniform dblRayGenData record) { SolidAccelerationStructure world = record.ray[rayID].volume_accel_solid; // Pass the ray's origin and direction to the payload - payload.primitive_id = xdg::ID_NONE; // Initialize to ID_NONE to indicate the point is outside all elements + payload.primitive_id = ID_NONE; // Initialize to -1 to indicate the point is TracePoint(world, RAY_FLAG_NONE, 0xff, xdg::RT_VOLUME_RAY_INDEX, xdg::RT_VOLUME_MISS_INDEX, pointDesc, payload); @@ -41,7 +41,7 @@ void find_element(uniform dblRayGenData record) { [shader("miss")] void tet_miss(inout SolidRayFirePayload payload) { - payload.primitive_id = xdg::ID_NONE; // Mark the ray as outside all tetrahedra + payload.primitive_id = ID_NONE; // Mark the ray as outside all tetrahedra } [shader("closesthit")] diff --git a/src/gprt/triangle_rt_shaders.slang b/src/gprt/triangle_rt_shaders.slang index 600f6411..f8e63ab3 100644 --- a/src/gprt/triangle_rt_shaders.slang +++ b/src/gprt/triangle_rt_shaders.slang @@ -32,8 +32,8 @@ void ray_fire_hit(uniform DPTriangleGeomData record, inout SurfaceRayFirePayload void ray_fire_miss(inout SurfaceRayFirePayload payload) { // Set the miss payload to default values payload.distance = -1.0f; - payload.surf_id = xdg::ID_NONE; - payload.primitive_id = xdg::ID_NONE; + payload.surf_id = ID_NONE; + payload.primitive_id = ID_NONE; payload.piv = xdg::PointInVolume::UNSET; } @@ -55,7 +55,7 @@ void ray_fire(uniform dblRayGenData record) { // Pass the ray's origin and direction to the payload payload.distance = -1.0f; - payload.surf_id = xdg::ID_NONE; + payload.surf_id = ID_NONE; payload.tlas = world; TraceRay(world, RAY_FLAG_NONE, 0xff, xdg::RT_SURFACE_RAY_INDEX, xdg::RT_SURFACE_MISS_INDEX, rayDesc, payload); @@ -81,7 +81,7 @@ void point_in_volume(uniform dblRayGenData record) { SurfaceAccelerationStructure world = record.ray[rayID].volume_accel_surf; // Pass the ray's origin and direction to the payload - payload.surf_id = xdg::ID_NONE; + payload.surf_id = ID_NONE; payload.tlas = world; payload.piv = xdg::PointInVolume::OUTSIDE; // Initialize point in volume check result to outside (0) diff --git a/vendor/GPRT b/vendor/GPRT index 405d9ee9..f1e95e41 160000 --- a/vendor/GPRT +++ b/vendor/GPRT @@ -1 +1 @@ -Subproject commit 405d9ee9f5ee8e1a0455f776f9e2c3adffb64160 +Subproject commit f1e95e4188cde591547d6b4a33a70bf2afaeec59