diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 393e6ad6..c5fe1de5 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -13,9 +13,9 @@ jobs: steps: - name: Checkout - uses: actions/checkout@v2 + uses: actions/checkout@v4 with: - submodules: true + submodules: recursive - name: Apt dependencies shell: bash @@ -27,7 +27,31 @@ jobs: libnetcdf-dev \ libhdf5-dev \ libeigen3-dev \ - cmake + cmake + - name: Update CMake + uses: ssrobins/install-cmake@v1 + with: + version: 3.24.3 + cache: true + + - name: GPRT Apt Dependencies + shell: bash + run: | + sudo apt install -y \ + libpthread-stubs0-dev \ + libx11-xcb-dev \ + xorg-dev \ + libxinerama-dev \ + libglu1-mesa-dev \ + freeglut3-dev \ + mesa-common-dev \ + libglfw3 + + - name: Prepare Vulkan SDK + uses: humbletim/install-vulkan-sdk@v1.2 + with: + version: 1.4.309.0 + cache: true - name: MOAB Clone shell: bash @@ -97,9 +121,10 @@ jobs: - name: Build shell: bash run: | + export LD_LIBRARY_PATH=$VULKAN_SDK/lib:$LD_LIBRARY_PATH mkdir build cd build - cmake .. -DCMAKE_PREFIX_PATH="$HOME/LIBMESH;$HOME/MOAB" -DCMAKE_INSTALL_PREFIX=$HOME/opt -DXDG_ENABLE_MOAB=ON -DXDG_ENABLE_LIBMESH=ON + cmake .. -DCMAKE_PREFIX_PATH="$HOME/LIBMESH;$HOME/MOAB" -DCMAKE_INSTALL_PREFIX=$HOME/opt -DXDG_ENABLE_MOAB=ON -DXDG_ENABLE_LIBMESH=ON -DXDG_ENABLE_GPRT=ON make -j4 install - name: Test diff --git a/.gitmodules b/.gitmodules index 771eab47..1455f373 100644 --- a/.gitmodules +++ b/.gitmodules @@ -17,3 +17,6 @@ [submodule "vendor/linalg"] path = vendor/linalg url = git@github.com:sgorsten/linalg.git +[submodule "vendor/GPRT"] + path = vendor/GPRT + url = https://github.com/gprt-org/GPRT.git diff --git a/CMakeLists.txt b/CMakeLists.txt index b9df5518..37b97734 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -5,6 +5,7 @@ option(XDG_ENABLE_MOAB "Enable support for the MOAB mesh library" ON option(XDG_ENABLE_MFEM "Enable support for the MFEM mesh library" OFF) option(XDG_ENABLE_LIBMESH "Enable support for the libMesh mesh library" OFF) option(XDG_ENABLE_EMBREE "Enable support for the Embree ray tracing library" ON) +option(XDG_ENABLE_GPRT "Enable support for the GPRT ray tracing library" OFF) option(XDG_BUILD_TESTS "Enable C++ unit testing" ON) option(XDG_BUILD_TOOLS "Enable tools and miniapps" ON) @@ -93,6 +94,21 @@ add_subdirectory(vendor/argparse) set(INDICATORS_INSTALL ON CACHE STRING "Ensure indicators targets are provided") add_subdirectory(vendor/indicators) +# Ensure at least one ray tracing backend is enabled +if (NOT XDG_ENABLE_EMBREE AND NOT XDG_ENABLE_GPRT) + message(FATAL_ERROR + "No ray tracing backend enabled. Enable at least one of:\n" + " -DXDG_ENABLE_EMBREE=ON\n" + " -DXDG_ENABLE_GPRT=ON") +endif() + + +# GPRT +set(GPRT_BUILD_SHARED ON CACHE BOOL "Build GPRT as a shared library" FORCE) +if (XDG_ENABLE_GPRT) + add_subdirectory(vendor/GPRT) +endif() + list(APPEND xdg_sources src/geometry/measure.cpp src/geometry/plucker.cpp @@ -114,6 +130,16 @@ src/embree/ray_tracer.cpp ) endif() +if (XDG_ENABLE_GPRT) +list(APPEND xdg_sources +src/gprt/ray_tracer.cpp +) +list(APPEND xdg_device_codes +dbl_deviceCode +) + +endif() + if (XDG_ENABLE_LIBMESH) list(APPEND xdg_sources src/libmesh/mesh_manager.cpp @@ -202,6 +228,11 @@ if (XDG_ENABLE_EMBREE) target_compile_definitions(xdg PUBLIC XDG_ENABLE_EMBREE) endif() +if (XDG_ENABLE_GPRT) + target_compile_definitions(xdg PUBLIC XDG_ENABLE_GPRT) + target_link_options(xdg INTERFACE -Wl,--unresolved-symbols=ignore-in-shared-libs) +endif() + # ========================== # Link ray tracing libraries # ========================== @@ -210,6 +241,22 @@ if (XDG_ENABLE_EMBREE) target_link_libraries(xdg embree fmt::fmt) endif() +if (XDG_ENABLE_GPRT) + # Compile device code for each slang file (currently only dbl_deviceCode.slang) + foreach(device_code ${xdg_device_codes}) + embed_devicecode( + OUTPUT_TARGET + ${device_code} + HEADERS + ${CMAKE_CURRENT_SOURCE_DIR}/include/xdg/gprt/shared_structs.h + SOURCES + ${CMAKE_CURRENT_SOURCE_DIR}/src/gprt/${device_code}.slang + ) + target_link_libraries(xdg ${device_code}) + endforeach() +target_link_libraries(xdg $) +endif() + # =================== # Link mesh libraries # =================== @@ -235,6 +282,24 @@ install(TARGETS xdg INCLUDES DESTINATION include ) +# Install GPRT Targets from shared library +if (TARGET gprt) + install(TARGETS gprt + RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR} + LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR} + ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR} + ) +endif() + +# Install targets for each slang file compiled (currently only dbl_deviceCode) +foreach(device_code ${xdg_device_codes}) + install(TARGETS ${device_code} + EXPORT xdg-targets + LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR} + ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR} + ) +endforeach() + install(EXPORT xdg-targets FILE XDGTargets.cmake NAMESPACE xdg:: diff --git a/include/xdg/constants.h b/include/xdg/constants.h index e1e3245d..b31230d6 100644 --- a/include/xdg/constants.h +++ b/include/xdg/constants.h @@ -6,6 +6,10 @@ #include #include +#include "fmt/format.h" + +#include "xdg/shared_enums.h" + namespace xdg { constexpr double INFTY {std::numeric_limits::max()}; @@ -63,7 +67,7 @@ static const std::map MESH_LIB_TO_STR = static const std::map RT_LIB_TO_STR = { {RTLibrary::EMBREE, "EMBREE"}, - {RTLibrary::GPRT, "GPRT"}, + {RTLibrary::GPRT, "GPRT"} }; // Mesh identifer type @@ -111,9 +115,6 @@ static Property VOID_MATERIAL {PropertyType::MATERIAL, "void"}; // Enumerator for different ray fire types enum class RayFireType { VOLUME, POINT_CONTAINMENT, ACCUMULATE_HITS, FIND_VOLUME }; -// -enum class HitOrientation { ANY, EXITING, ENTERING }; - // Enumerator for different element types (maybe we want more here?) enum class SurfaceElementType { TRI = 0, @@ -127,4 +128,22 @@ enum class VolumeElementType { } // namespace xdg +namespace fmt { + template <> +struct formatter : fmt::formatter { + auto format(xdg::RTLibrary lib, fmt::format_context& ctx) { + return fmt::formatter::format(xdg::RT_LIB_TO_STR.at(lib), ctx); + } +}; + +template <> +struct formatter : fmt::formatter { + auto format(xdg::MeshLibrary lib, fmt::format_context& ctx) { + return fmt::formatter::format(xdg::MESH_LIB_TO_STR.at(lib), ctx); + } +}; + + +} + #endif // include guard \ No newline at end of file diff --git a/include/xdg/embree/ray_tracer.h b/include/xdg/embree/ray_tracer.h index a35afdf3..4ec064b5 100644 --- a/include/xdg/embree/ray_tracer.h +++ b/include/xdg/embree/ray_tracer.h @@ -21,6 +21,8 @@ class EmbreeRayTracer : public RayTracer { public: EmbreeRayTracer(); ~EmbreeRayTracer(); + RTLibrary library() const override { return RTLibrary::EMBREE; } + void init() override; RTCScene create_embree_scene(); @@ -38,6 +40,7 @@ class EmbreeRayTracer : public RayTracer { MeshID find_element(TreeID tree, const Position& point) const override; + // Query Methods bool point_in_volume(TreeID scene, const Position& point, @@ -88,6 +91,7 @@ class EmbreeRayTracer : public RayTracer { MeshID surface, RTCScene& volume_scene, int& storage_offset); + // Global Tree IDs RTCScene global_surface_scene_ {nullptr}; RTCScene global_element_scene_ {nullptr}; diff --git a/include/xdg/gprt/ray_tracer.h b/include/xdg/gprt/ray_tracer.h new file mode 100644 index 00000000..0805af27 --- /dev/null +++ b/include/xdg/gprt/ray_tracer.h @@ -0,0 +1,132 @@ +#ifndef _XDG_GPRT_BASE_RAY_TRACING_INTERFACE_H +#define _XDG_GPRT_BASE_RAY_TRACING_INTERFACE_H + +#include +#include +#include + +#include "xdg/constants.h" +#include "xdg/mesh_manager_interface.h" +#include "xdg/primitive_ref.h" +#include "xdg/geometry_data.h" +#include "xdg/ray_tracing_interface.h" +#include "xdg/ray.h" +#include "xdg/error.h" +#include "gprt/gprt.h" +#include "shared_structs.h" + +extern GPRTProgram dbl_deviceCode; +namespace xdg { + +class GPRTRayTracer : public RayTracer { + public: + GPRTRayTracer(); + ~GPRTRayTracer(); + RTLibrary library() const override { return RTLibrary::GPRT; } + + void set_geom_data(const std::shared_ptr mesh_manager); + void init() override; + + // 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(TreeID tree, const Position& point) const override { + fatal_error("Element trees not currently supported with GPRT ray tracer"); + return ID_NONE; + }; + + std::pair + register_volume(const std::shared_ptr& mesh_manager, + MeshID volume) override; + + TreeID create_surface_tree(const std::shared_ptr& mesh_manager, + MeshID volume) override; + + TreeID create_element_tree(const std::shared_ptr& mesh_manager, + MeshID volume) override; + + void create_global_surface_tree() override; + + void create_global_element_tree() override + { + warning("Global element trees not currently supported with GPRT ray tracer"); + return; + }; + + bool point_in_volume(TreeID scene, + const Position& point, + const Direction* direction = nullptr, + const std::vector* exclude_primitives = nullptr) const override; + + std::pair ray_fire(TreeID scene, + const Position& origin, + const Direction& direction, + const double dist_limit = INFTY, + HitOrientation orientation = HitOrientation::EXITING, + std::vector* const exclude_primitives = nullptr) override; + + void closest(TreeID scene, + const Position& origin, + double& dist, + MeshID& triangle) override {}; + void closest(TreeID scene, + const Position& origin, + double& dist) override {}; + + bool occluded(TreeID scene, + const Position& origin, + const Direction& direction, + double& dist) const override { + fatal_error("Occlusion queries are not currently supported with GPRT ray tracer"); + return false; + } + + private: + // GPRT objects + GPRTContext context_; + GPRTProgram deviceCode_; // device code for float precision shaders + GPRTModule module_; // device code module for single precision shaders + GPRTAccel world_; + GPRTBuildParams buildParams_; // rayGenProgram_; + GPRTRayGenOf rayGenPointInVolProgram_; + GPRTMissOf missProgram_; + GPRTComputeOf aabbPopulationProgram_; // rayInputBuffer_; // rayOutputBuffer_; // excludePrimitivesBuffer_; // globalBlasInstances_; // trianglesGeomType_; //> surface_to_geometry_map_; // surface_volume_tree_to_accel_map; // Map from XDG::TreeID to GPRTAccel for volume TLAS + std::vector blas_handles_; // Store BLAS handles so that they can be explicitly referenced in destructor + + // Global Tree IDs + GPRTAccel global_surface_accel_ {nullptr}; + GPRTAccel global_element_accel_ {nullptr}; + + }; + +} // namespace xdg + +#endif // include guard \ No newline at end of file diff --git a/include/xdg/gprt/shared_structs.h b/include/xdg/gprt/shared_structs.h new file mode 100644 index 00000000..09192988 --- /dev/null +++ b/include/xdg/gprt/shared_structs.h @@ -0,0 +1,60 @@ +#include "gprt.h" +#include "../shared_enums.h" + +struct GPRTPrimitiveRef +{ + int id; // ID of the primitive + int sense; +}; + +struct dblRayInput +{ + double3 origin; + double3 direction; + double tMin; // Minimum distance for ray intersection + double tMax; // Maximum distance for ray intersection + int32_t* exclude_primitives; // Optional for excluding primitives + 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 +}; + +struct dblRayOutput +{ + double distance; + int surf_id; + int primitive_id; + xdg::PointInVolume piv; // Point in volume check result (0 for outside, 1 for inside) +}; + +/* variables for double precision triangle mesh geometry */ +struct DPTriangleGeomData { + double3 *vertex; // vertex buffer + float3 *aabbs; // AABB buffer + uint3 *index; // index buffer + double3 *normals; // normals buffer + int surf_id; + int2 vols; + int forward_vol; + int reverse_vol; + dblRayInput *rayIn; // double precision rays + xdg::HitOrientation hitOrientation; + int forward_tree; // TreeID of the forward volume + int reverse_tree; // TreeID of the reverse volume + GPRTPrimitiveRef* primitive_refs; + int num_faces; // Number of faces in the geometry +}; + +struct dblRayGenData { + dblRayInput *ray; + dblRayOutput *out; +}; + +/* A small structure of constants that can change every frame without rebuilding the + shader binding table. (must be 128 bytes or less) */ + +struct dblRayFirePushConstants { + double tMax; + double tMin; +}; diff --git a/include/xdg/mesh_managers.h b/include/xdg/mesh_managers.h new file mode 100644 index 00000000..fe6912ce --- /dev/null +++ b/include/xdg/mesh_managers.h @@ -0,0 +1,8 @@ +// mesh manager concrete implementations +#ifdef XDG_ENABLE_MOAB +#include "xdg/moab/mesh_manager.h" +#endif + +#ifdef XDG_ENABLE_LIBMESH +#include "xdg/libmesh/mesh_manager.h" +#endif diff --git a/include/xdg/ray_tracers.h b/include/xdg/ray_tracers.h new file mode 100644 index 00000000..b4816571 --- /dev/null +++ b/include/xdg/ray_tracers.h @@ -0,0 +1,8 @@ +// ray tracing interface concrete implementations +#ifdef XDG_ENABLE_EMBREE +#include "xdg/embree/ray_tracer.h" +#endif + +#ifdef XDG_ENABLE_GPRT +#include "xdg/gprt/ray_tracer.h" +#endif \ No newline at end of file diff --git a/include/xdg/ray_tracing_interface.h b/include/xdg/ray_tracing_interface.h index 4b50c922..3ea7470d 100644 --- a/include/xdg/ray_tracing_interface.h +++ b/include/xdg/ray_tracing_interface.h @@ -121,6 +121,9 @@ class RayTracer { const Direction& direction, double& dist) const = 0; + virtual RTLibrary library() const = 0; + + // Generic Accessors int num_registered_trees() const { return surface_trees_.size() + element_trees_.size(); }; int num_registered_surface_trees() const { return surface_trees_.size(); }; diff --git a/include/xdg/shared_enums.h b/include/xdg/shared_enums.h new file mode 100644 index 00000000..f6198f60 --- /dev/null +++ b/include/xdg/shared_enums.h @@ -0,0 +1,19 @@ +#ifndef XDG_SHARED_ENUMS_H +#define XDG_SHARED_ENUMS_H + +namespace xdg { + + enum PointInVolume : int { + 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/include/xdg/xdg.h b/include/xdg/xdg.h index 2d7f838c..fe7c3733 100644 --- a/include/xdg/xdg.h +++ b/include/xdg/xdg.h @@ -6,7 +6,7 @@ #include "xdg/mesh_manager_interface.h" #include "xdg/ray_tracing_interface.h" -#include "xdg/embree/ray_tracer.h" + namespace xdg { @@ -16,19 +16,7 @@ class XDG { // Constructors XDG() = default; - XDG(std::shared_ptr mesh_manager, RTLibrary ray_tracing_lib = RTLibrary::EMBREE) : mesh_manager_(mesh_manager) - { - // construct internal raytracer for XDG - switch (ray_tracing_lib) - { - case RTLibrary::EMBREE: - set_ray_tracing_interface(std::make_shared()); - break; - case RTLibrary::GPRT: - fatal_error("This backend is not yet implemented"); - break; - } - } + XDG(std::shared_ptr mesh_manager, RTLibrary ray_tracing_lib = RTLibrary::EMBREE); // factory method that allows for specification of a backend mesh library and ray tracer. Default to MOAB + EMBREE static std::shared_ptr create(MeshLibrary mesh_lib = MeshLibrary::MOAB, RTLibrary ray_tracing_lib = RTLibrary::EMBREE); @@ -140,7 +128,6 @@ Direction surface_normal(MeshID surface, std::unordered_map volume_to_surface_tree_map_; // surface_to_tree_map_; // volume_to_point_location_tree_map_; // surface_to_geometry_map_; // vertices, return EXIT_EARLY; } - // Determine the value of the second Plucker coordinate from edge 2 + // Determine the value of the third Plucker coordinate from edge 2 double plucker_coord2 = plucker_edge_test(vertices[2], vertices[0], raya, rayb); @@ -96,8 +96,8 @@ bool plucker_ray_tri_intersect(const std::array vertices, // get the distance to intersection const double inverse_sum = 1.0 / (plucker_coord0 + plucker_coord1 + plucker_coord2); - // TODO: replace assert with warning assert(0.0 != inverse_sum); + const Position intersection(plucker_coord0 * inverse_sum * vertices[2] + plucker_coord1 * inverse_sum * vertices[0] + plucker_coord2 * inverse_sum * vertices[1]); @@ -111,17 +111,19 @@ bool plucker_ray_tri_intersect(const std::array vertices, max_abs_dir = fabs(direction[i]); } } + dist_out = (intersection[idx] - origin[idx]) / direction[idx]; // is the intersection within distance limits? - if( ( nonneg_ray_len && nonneg_ray_len < dist_out ) || // intersection is beyond positive limit - ( neg_ray_len && *neg_ray_len >= dist_out ) || // intersection is behind negative limit - ( !neg_ray_len && 0 > dist_out ) ) - { // Unless a neg_ray_len is used, don't return negative distances - return EXIT_EARLY; + if ((nonneg_ray_len && nonneg_ray_len < dist_out) || // intersection is beyond positive limit + (neg_ray_len && *neg_ray_len >= dist_out) || // intersection is behind negative limit + (!neg_ray_len && 0 > dist_out)) // unless neg_ray_len used, don't allow negative distances + { + return EXIT_EARLY; } return true; } + } // namespace xdg \ No newline at end of file diff --git a/src/gprt/dbl_deviceCode.slang b/src/gprt/dbl_deviceCode.slang new file mode 100644 index 00000000..2f069704 --- /dev/null +++ b/src/gprt/dbl_deviceCode.slang @@ -0,0 +1,316 @@ +#include "../../include/xdg/gprt/shared_structs.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 rayOutput 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.out[rayID].distance = payload.distance; + record.out[rayID].surf_id = payload.surf_id; + record.out[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.out.surf_id = payload.surf_id; + record.out[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*/ +[shader("compute")] +[numthreads(1, 1, 1)] +void +populate_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; + uint flags = RayFlags(); + + if (rayID >= nRays) { + return; + } + + bool useOrientation = false; + int orientation = 0; + if ((flags & RAY_FLAG_CULL_BACK_FACING_TRIANGLES) != 0) { + orientation = -1; + useOrientation = true; + } + else if ((flags & RAY_FLAG_CULL_FRONT_FACING_TRIANGLES) != 0) { + orientation = 1; + useOrientation = true; + } + + int3 indices = record.index[primID]; + double3 v0 = record.vertex[indices[0]]; + double3 v1 = record.vertex[indices[1]]; + double3 v2 = record.vertex[indices[2]]; + + double3 origin = record.rayIn[rayID].origin; + double3 direction = record.rayIn[rayID].direction; + + // double tMin = record.rayIn[rayID].tMin; + double tMin = record.rayIn[rayID].tMin; + double tMax = record.rayIn[rayID].tMax; + + const double3 raya = direction; + const double3 rayb = cross(direction, origin); + + double plucker_coord0 = plucker_edge_test(v0, v1, raya, rayb); + if (useOrientation && orientation * plucker_coord0 > 0) { + return; + } + + double plucker_coord1 = plucker_edge_test(v1, v2, raya, rayb); + if (useOrientation && orientation * plucker_coord1 > 0) { + return; + } + else if ((0.0 < plucker_coord0 && 0.0 > plucker_coord1) || (0.0 > plucker_coord0 && 0.0 < plucker_coord1)) { + return; + } + + double plucker_coord2 = plucker_edge_test(v2, v0, raya, rayb); + if (useOrientation && orientation * plucker_coord2 > 0) { + return; + } + else if ((0.0 < plucker_coord1 && 0.0 > plucker_coord2) || (0.0 > plucker_coord1 && 0.0 < plucker_coord2) || + (0.0 < plucker_coord0 && 0.0 > plucker_coord2) || (0.0 > plucker_coord0 && 0.0 < plucker_coord2)) { + return; + } + + if (0.0 == plucker_coord0 && 0.0 == plucker_coord1 && 0.0 == plucker_coord2) { + return; + } + + const double inverse_sum = 1.0 / (plucker_coord0 + plucker_coord1 + plucker_coord2); + const double3 intersection = double3(plucker_coord0 * inverse_sum * v2 + + plucker_coord1 * inverse_sum * v0 + + plucker_coord2 * inverse_sum * v1); + + int idx = 0; + double max_abs_dir = 0; + for(uint i = 0; i < 3; ++i) { + if(abs(direction[i]) > max_abs_dir) { + idx = i; + max_abs_dir = abs(direction[i]); + } + } + const double dist = (intersection[idx] - origin[idx]) / direction[idx]; + + double t = dist; + double u = plucker_coord2 * inverse_sum; + double v = plucker_coord0 * inverse_sum; + + + if(u < 0.0 || v < 0.0 || (u + v) > 1.0) { + t = -1.0; + } + if (t > tMax) { + return; + } + if (t < tMin) { + return; + } + + DPAttribute attr; + attr.f64t = t; + + float f32t = float(t); + if (double(f32t) < t) f32t = next_after(f32t); + + double3 norm = record.normals[primID]; // recover double precision normal. TODO - Should we calculate from vertices instead? + + // sense adjustment of normal + if (record.rayIn[rayID].volume_tree == record.reverse_tree) + { + norm = -norm; + } + + double norm_dot_dir = dot(norm, direction); + uint hit_kind = norm_dot_dir < 0 ? HIT_KIND_TRIANGLE_FRONT_FACE + : HIT_KIND_TRIANGLE_BACK_FACE; + + xdg::HitOrientation hitOrientation = record.rayIn[rayID].hitOrientation; + + if (orientation_cull(direction, norm, hitOrientation)) + { + return; + } + + for (int i = 0; i < record.rayIn[rayID].exclude_count; ++i) + { + if (record.rayIn[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 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 +} + +// Plucker coordinate +double plucker_edge_test(in double3 vertexa, in double3 vertexb, in double3 ray, in double3 ray_normal) +{ + double pip; + const double near_zero = 10 * DBL_EPSILON; + + if (first(vertexa, vertexb)) + { + double3 edge = vertexb - vertexa; + double3 edge_normal = cross(edge, vertexa); + pip = dot(ray, edge_normal) + dot(ray_normal, edge); + } + else + { + double3 edge = vertexa - vertexb; + double3 edge_normal = cross(edge, vertexb); + pip = dot(ray, edge_normal) + 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 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_); +} \ No newline at end of file diff --git a/src/gprt/ray_tracer.cpp b/src/gprt/ray_tracer.cpp new file mode 100644 index 00000000..c5778f97 --- /dev/null +++ b/src/gprt/ray_tracer.cpp @@ -0,0 +1,337 @@ +#include "xdg/gprt/ray_tracer.h" +#include "gprt/gprt.h" + +namespace xdg { + +GPRTRayTracer::GPRTRayTracer() +{ + gprtRequestRayTypeCount(numRayTypes_); // Set the number of shaders which can be set to the same geometry + context_ = gprtContextCreate(); + module_ = gprtModuleCreate(context_, dbl_deviceCode); + + numRays = 1; // Set the number of rays to be cast + rayInputBuffer_ = gprtDeviceBufferCreate(context_, numRays); + rayOutputBuffer_ = gprtDeviceBufferCreate(context_, numRays); + excludePrimitivesBuffer_ = gprtDeviceBufferCreate(context_); // initialise buffer of size 1 + + setup_shaders(); + + // Bind the buffers to the RayGenData structure + dblRayGenData* rayGenData = gprtRayGenGetParameters(rayGenProgram_); + rayGenData->ray = gprtBufferGetDevicePointer(rayInputBuffer_); + rayGenData->out = gprtBufferGetDevicePointer(rayOutputBuffer_); + + // Bind the buffers to the RayGenData structure + dblRayGenData* rayGenPIVData = gprtRayGenGetParameters(rayGenPointInVolProgram_); + rayGenPIVData->ray = gprtBufferGetDevicePointer(rayInputBuffer_); + rayGenPIVData->out = gprtBufferGetDevicePointer(rayOutputBuffer_); + + // Set up build parameters for acceleration structures + buildParams_.buildMode = GPRT_BUILD_MODE_FAST_BUILD_NO_UPDATE; +} + +GPRTRayTracer::~GPRTRayTracer() +{ + // Ensure all GPU operations are complete before destroying resources + gprtGraphicsSynchronize(context_); + gprtComputeSynchronize(context_); + + + // Destroy TLAS structures + for (const auto& [tree, accel] : surface_volume_tree_to_accel_map) { + gprtAccelDestroy(accel); + } + + // Destroy BLAS structures + for (const auto& blas : blas_handles_) { + gprtAccelDestroy(blas); + } + + // Destroy Geoms and Types + for (const auto& [surf, geom] : surface_to_geometry_map_) { + gprtGeomDestroy(geom); + } + gprtGeomTypeDestroy(trianglesGeomType_); + + // Destroy Buffers + gprtBufferDestroy(rayInputBuffer_); + gprtBufferDestroy(rayOutputBuffer_); + gprtBufferDestroy(excludePrimitivesBuffer_); + + // Destroy module and context + gprtModuleDestroy(module_); + gprtContextDestroy(context_); +} + +void GPRTRayTracer::setup_shaders() +{ + // Set up ray generation and miss programs + rayGenProgram_ = gprtRayGenCreate(context_, module_, "ray_fire"); + rayGenPointInVolProgram_ = gprtRayGenCreate(context_, module_, "point_in_volume"); + missProgram_ = gprtMissCreate(context_, module_, "ray_fire_miss"); + aabbPopulationProgram_ = gprtComputeCreate(context_, module_, "populate_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 +} + +void GPRTRayTracer::init() +{ + // Build the shader binding table (SBT) after all shader programs and acceleration structures are set up + gprtBuildShaderBindingTable(context_, GPRT_SBT_ALL); + // Note that should we need to update any shaders or acceleration structures, we must rebuild the SBT +} + +std::pair +GPRTRayTracer::register_volume(const std::shared_ptr& mesh_manager, MeshID volume_id) +{ + // set up ray tracing tree for boundary faces of the volume + TreeID faces_tree = create_surface_tree(mesh_manager, volume_id); + // set up point location tree for any volumetric elements. TODO - currently not supported with GPRT + TreeID element_tree = create_element_tree(mesh_manager, volume_id); + return {faces_tree, element_tree}; // return TREE_NONE for element tree until implmemented +} + +SurfaceTreeID +GPRTRayTracer::create_surface_tree(const std::shared_ptr& mesh_manager, MeshID volume_id) +{ + SurfaceTreeID tree = next_surface_tree_id(); + surface_trees_.push_back(tree); + auto volume_surfaces = mesh_manager->get_volume_surfaces(volume_id); + 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 + + // Get storage for vertices + auto [vertices, indices] = mesh_manager->get_surface_mesh(surf); + 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 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]); + } + + // 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 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); + } + + auto vertex_buffer = gprtDeviceBufferCreate(context_, dbl3Vertices.size(), dbl3Vertices.data()); + auto aabb_buffer = gprtDeviceBufferCreate(context_, 2*num_faces, 0); // AABBs for each triangle + gprtAABBsSetPositions(triangleGeom, aabb_buffer, num_faces, 2*sizeof(float3), 0); + auto connectivity_buffer = gprtDeviceBufferCreate(context_, ui3Indices.size(), ui3Indices.data()); + auto normal_buffer = gprtDeviceBufferCreate(context_, num_faces, normals.data()); + auto primitive_refs_buffer = gprtDeviceBufferCreate(context_, num_faces, primitive_refs.data()); // Buffer for primitive sense + + geom_data->vertex = gprtBufferGetDevicePointer(vertex_buffer); + geom_data->index = gprtBufferGetDevicePointer(connectivity_buffer); + geom_data->aabbs = gprtBufferGetDevicePointer(aabb_buffer); + geom_data->rayIn = gprtBufferGetDevicePointer(rayInputBuffer_); + geom_data->surf_id = surf; + geom_data->normals = gprtBufferGetDevicePointer(normal_buffer); + 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); + gprtComputeSynchronize(context_); // Ensure all GPU operations are complete before accessing results + + GPRTAccel blas = gprtAABBAccelCreate(context_, triangleGeom, buildParams_.buildMode); + + gprtAccelBuild(context_, blas, buildParams_); + + gprt::Instance 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 + surface_to_geometry_map_[surf] = triangleGeom; + + geom_data = gprtGeomGetParameters(triangleGeom); + instance = gprtAccelGetInstance(blas); + instance.mask = 0xff; + surfaceBlasInstances.push_back(instance); + globalBlasInstances_.push_back(instance); + + // Always update per-volume info + auto [forward_parent, reverse_parent] = mesh_manager->get_parent_volumes(surf); + if (volume_id == forward_parent) { + geom_data->forward_vol = forward_parent; + geom_data->forward_tree = tree; + } else if (volume_id == reverse_parent) { + geom_data->reverse_vol = reverse_parent; + geom_data->reverse_tree = tree; + } else { + fatal_error("Volume {} is not a parent of surface {}", volume_id, surf); + } + } + + // Create a TLAS (Top-Level Acceleration Structure) for all BLAS instances in this volume + auto instanceBuffer = gprtDeviceBufferCreate(context_, surfaceBlasInstances.size(), surfaceBlasInstances.data()); + GPRTAccel volume_tlas = gprtInstanceAccelCreate(context_, surfaceBlasInstances.size(), instanceBuffer); + gprtAccelBuild(context_, volume_tlas, buildParams_); + surface_volume_tree_to_accel_map[tree] = volume_tlas; + + return tree; +} + +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; +}; + +bool GPRTRayTracer::point_in_volume(SurfaceTreeID tree, + const Position& point, + const Direction* direction, + const std::vector* exclude_primitives) const +{ + GPRTAccel volume = surface_volume_tree_to_accel_map.at(tree); + dblRayGenData* rayGenPIVData = gprtRayGenGetParameters(rayGenPointInVolProgram_); + + // Use provided direction or if Direction == nulptr use default direction + Direction directionUsed = (direction != nullptr) ? Direction{direction->x, direction->y, direction->z} + : Direction{1. / std::sqrt(2.0), 1. / std::sqrt(2.0), 0.0}; + + gprtBufferMap(rayInputBuffer_); // Update the ray input buffer + dblRayInput* rayInput = gprtBufferGetHostPointer(rayInputBuffer_); + rayInput[0].volume_accel = gprtAccelGetDeviceAddress(volume); + rayInput[0].origin = {point.x, point.y, point.z}; + rayInput[0].direction = {directionUsed.x, directionUsed.y, directionUsed.z}; + rayInput[0].tMax = INFTY; // Set a large distance limit + rayInput[0].tMin = 0.0; + rayInput[0].volume_tree = tree; // Set the TreeID of the volume being queried + rayInput[0].hitOrientation = HitOrientation::ANY; // No orientation culling for point-in-volume check + + if (exclude_primitives) { + if (!exclude_primitives->empty()) gprtBufferResize(context_, excludePrimitivesBuffer_, exclude_primitives->size(), false); + gprtBufferMap(excludePrimitivesBuffer_); + std::copy(exclude_primitives->begin(), exclude_primitives->end(), gprtBufferGetHostPointer(excludePrimitivesBuffer_)); + gprtBufferUnmap(excludePrimitivesBuffer_); + + rayInput[0].exclude_primitives = gprtBufferGetDevicePointer(excludePrimitivesBuffer_); + rayInput[0].exclude_count = exclude_primitives->size(); + } + else { + // If no primitives are excluded, set the pointer to null and count to 0 + rayInput[0].exclude_primitives = nullptr; + rayInput[0].exclude_count = 0; + } + gprtBufferUnmap(rayInputBuffer_); // required to sync buffer back on GPU? + + gprtRayGenLaunch1D(context_, rayGenPointInVolProgram_, 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 output from the ray output buffer + gprtBufferMap(rayOutputBuffer_); + dblRayOutput* rayOutput = gprtBufferGetHostPointer(rayOutputBuffer_); + auto surface = rayOutput[0].surf_id; + auto piv = rayOutput[0].piv; // Point in volume check result + gprtBufferUnmap(rayOutputBuffer_); // 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 ray hit nothing, the point is outside volume + if (surface == ID_NONE) return false; + + return piv; +} + + +// This will launch the rays and run our shaders in the ray tracing pipeline +// miss shader returns dist = 0.0 and elementID = -1 +// closest hit shader returns dist = distance to hit and elementID = triangle ID +std::pair GPRTRayTracer::ray_fire(SurfaceTreeID tree, + const Position& origin, + const Direction& direction, + double dist_limit, + HitOrientation orientation, + std::vector* const exclude_primitives) +{ + GPRTAccel volume = surface_volume_tree_to_accel_map.at(tree); + dblRayGenData* rayGenData = gprtRayGenGetParameters(rayGenProgram_); + + gprtBufferMap(rayInputBuffer_); // Update the ray input buffer + dblRayInput* rayInput = gprtBufferGetHostPointer(rayInputBuffer_); + rayInput[0].volume_accel = gprtAccelGetDeviceAddress(volume); + rayInput[0].origin = {origin.x, origin.y, origin.z}; + rayInput[0].direction = {direction.x, direction.y, direction.z}; + rayInput[0].tMax = dist_limit; + rayInput[0].tMin = 0.0; + rayInput[0].hitOrientation = orientation; // Set orientation for the ray + rayInput[0].volume_tree = tree; // Set the TreeID of the volume being queried + + if (exclude_primitives) { + if (!exclude_primitives->empty()) gprtBufferResize(context_, excludePrimitivesBuffer_, exclude_primitives->size(), false); + gprtBufferMap(excludePrimitivesBuffer_); + std::copy(exclude_primitives->begin(), exclude_primitives->end(), gprtBufferGetHostPointer(excludePrimitivesBuffer_)); + gprtBufferUnmap(excludePrimitivesBuffer_); + + rayInput[0].exclude_primitives = gprtBufferGetDevicePointer(excludePrimitivesBuffer_); + rayInput[0].exclude_count = exclude_primitives->size(); + } + else { + // If no primitives are excluded, set the pointer to null and count to 0 + rayInput[0].exclude_primitives = nullptr; + rayInput[0].exclude_count = 0; + } + gprtBufferUnmap(rayInputBuffer_); // required to sync buffer back on GPU? + + gprtRayGenLaunch1D(context_, rayGenProgram_, 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 output from the ray output buffer + gprtBufferMap(rayOutputBuffer_); + dblRayOutput* rayOutput = gprtBufferGetHostPointer(rayOutputBuffer_); + auto distance = rayOutput[0].distance; + auto surface = rayOutput[0].surf_id; + auto primitive_id = rayOutput[0].primitive_id; + gprtBufferUnmap(rayOutputBuffer_); // 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 (surface == ID_NONE) + return {INFTY, ID_NONE}; + else + if (exclude_primitives) exclude_primitives->push_back(primitive_id); + return {distance, surface}; +} + +void GPRTRayTracer::create_global_surface_tree() +{ + // Create a TLAS (Top-Level Acceleration Structure) for all the volumes + auto globalBuffer = gprtDeviceBufferCreate(context_, globalBlasInstances_.size(), globalBlasInstances_.data()); + GPRTAccel global_accel = gprtInstanceAccelCreate(context_, globalBlasInstances_.size(), globalBuffer); + gprtAccelBuild(context_, global_accel, buildParams_); + + SurfaceTreeID tree = next_surface_tree_id(); + surface_trees_.push_back(tree); + surface_volume_tree_to_accel_map[tree] = global_accel; + global_surface_tree_ = tree; + global_surface_accel_ = global_accel; +} + +} // namespace xdg \ No newline at end of file diff --git a/src/mesh_manager_interface.cpp b/src/mesh_manager_interface.cpp index 4c762ee0..8a558354 100644 --- a/src/mesh_manager_interface.cpp +++ b/src/mesh_manager_interface.cpp @@ -263,5 +263,4 @@ MeshManager::get_parent_volumes(MeshID surface) const return this->surface_senses(surface); } - } // namespace xdg \ No newline at end of file diff --git a/src/xdg.cpp b/src/xdg.cpp index 63542bdd..f215917c 100644 --- a/src/xdg.cpp +++ b/src/xdg.cpp @@ -3,22 +3,37 @@ #include "xdg/xdg.h" #include "xdg/error.h" -#include "xdg/embree/ray_tracer.h" -// #include "xdg/gprt/ray_tracer.h" Not implemented yet +#include "xdg/constants.h" +#include "xdg/geometry/measure.h" -// mesh manager concrete implementations -#ifdef XDG_ENABLE_MOAB -#include "xdg/moab/mesh_manager.h" -#endif +#include "xdg/mesh_managers.h" -#ifdef XDG_ENABLE_LIBMESH -#include "xdg/libmesh/mesh_manager.h" -#endif +#include "xdg/ray_tracers.h" -#include "xdg/constants.h" -#include "xdg/geometry/measure.h" namespace xdg { +XDG::XDG(std::shared_ptr mesh_manager, RTLibrary ray_tracing_lib) + : mesh_manager_(mesh_manager) +{ + switch (ray_tracing_lib) { + case RTLibrary::EMBREE: + #ifdef XDG_ENABLE_EMBREE + set_ray_tracing_interface(std::make_shared()); + break; + #else + fatal_error("This build was not compiled with Embree support (XDG_ENABLE_EMBREE=OFF)."); + #endif + + case RTLibrary::GPRT: + #ifdef XDG_ENABLE_GPRT + set_ray_tracing_interface(std::make_shared()); + break; + #else + fatal_error("This build was not compiled with GPRT support (XDG_ENABLE_GPRT=OFF)."); + #endif + } +} + void XDG::prepare_raytracer() { for (auto volume : mesh_manager()->volumes()) { @@ -27,6 +42,8 @@ void XDG::prepare_raytracer() ray_tracing_interface()->create_global_element_tree(); ray_tracing_interface()->create_global_surface_tree(); + + ray_tracing_interface()->init(); // Initialize the ray tracer (e.g. build SBT for GPRT) } void XDG::prepare_volume_for_raytracing(MeshID volume) { @@ -65,7 +82,7 @@ std::shared_ptr XDG::create(MeshLibrary mesh_lib, RTLibrary ray_tracing_lib if (ray_tracing_lib == RTLibrary::EMBREE) return std::make_shared(); #endif #ifdef XDG_ENABLE_GPRT - // if (ray_tracing_lib == RTLibrary::GPRT) return std::make_shared(); + if (ray_tracing_lib == RTLibrary::GPRT) return std::make_shared(); #endif // If no supported ray tracing library throw an error diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 2f537135..3e8631c7 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -41,7 +41,9 @@ endif() foreach(test ${TEST_NAMES}) add_executable(${test} ${test}.cpp) target_link_libraries(${test} xdg Catch2::Catch2WithMain) - add_test(NAME ${test} COMMAND ${test} WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) + add_test(NAME ${test} + COMMAND ${test} --allow-running-no-tests + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) endforeach() set( diff --git a/tests/test_moab.cpp b/tests/test_moab.cpp index 921d4d43..4638a3f7 100644 --- a/tests/test_moab.cpp +++ b/tests/test_moab.cpp @@ -5,13 +5,14 @@ // testing includes #include #include +#include // xdg includes #include "xdg/error.h" #include "xdg/mesh_manager_interface.h" #include "xdg/moab/mesh_manager.h" #include "xdg/xdg.h" -#include "xdg/embree/ray_tracer.h" +#include "util.h" using namespace xdg; @@ -69,47 +70,56 @@ TEST_CASE("Test BVH Build") REQUIRE(mesh_manager->num_volumes() == 2); REQUIRE(mesh_manager->num_surfaces() == 6); - std::unique_ptr ray_tracing_interface = std::make_unique(); + // Generate one test run per enabled backend + auto rt_backend = GENERATE(RTLibrary::EMBREE, RTLibrary::GPRT); - for (auto volume : mesh_manager->volumes()) { - ray_tracing_interface->register_volume(mesh_manager, volume); - } - - REQUIRE(ray_tracing_interface->num_registered_trees() == 2); + // Actual testing logic + DYNAMIC_SECTION(fmt::format("Backend = {}", rt_backend)) { + check_ray_tracer_supported(rt_backend); // skip if backend not enabled at configuration time + auto rti = create_raytracer(rt_backend); + + for (const auto& volume : mesh_manager->volumes()) { + rti->register_volume(mesh_manager, volume); + } + REQUIRE(rti->num_registered_trees() == 2); + } } -TEST_CASE("Test Ray Fire MOAB") -{ - std::shared_ptr xdg = XDG::create(MeshLibrary::MOAB); - REQUIRE(xdg->mesh_manager()->mesh_library() == MeshLibrary::MOAB); - const auto& mesh_manager = xdg->mesh_manager(); - mesh_manager->load_file("cube.h5m"); - mesh_manager->init(); - xdg->prepare_raytracer(); +TEST_CASE("Test Ray Fire MOAB (all built backends)", "[ray_tracer][moab]") +{ + // Generate one test run per enabled backend + auto rt_backend = GENERATE(RTLibrary::EMBREE, RTLibrary::GPRT); - MeshID volume = mesh_manager->volumes()[0]; + DYNAMIC_SECTION(fmt::format("Backend = {}", rt_backend)) { + check_ray_tracer_supported(rt_backend); // skip if backend not enabled at configuration time + auto xdg = XDG::create(MeshLibrary::MOAB, rt_backend); + REQUIRE(xdg->mesh_manager()->mesh_library() == MeshLibrary::MOAB); - Position origin {0.0, 0.0, 0.0}; - Direction direction {1.0, 0.0, 0.0}; - std::pair intersection; + const auto& mm = xdg->mesh_manager(); + mm->load_file("cube.h5m"); + mm->init(); + xdg->prepare_raytracer(); - intersection = xdg->ray_fire(volume, origin, direction); + MeshID volume = mm->volumes()[0]; - // this cube is 10 cm on a side, so the ray should hit the surface at 5 cm - REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(5.0, 1e-6)); + Position origin {0.0, 0.0, 0.0}; + Direction dir {1.0, 0.0, 0.0}; - origin = {3.0, 0.0, 0.0}; - intersection = xdg->ray_fire(volume, origin, direction); - REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(2.0, 1e-6)); + auto hit = xdg->ray_fire(volume, origin, dir); + REQUIRE_THAT(hit.first, Catch::Matchers::WithinAbs(5.0, 1e-6)); - origin = {-10.0, 0.0, 0.0}; - intersection = xdg->ray_fire(volume, origin, direction); - REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(15.0, 1e-6)); + origin = {3.0, 0.0, 0.0}; + hit = xdg->ray_fire(volume, origin, dir); + REQUIRE_THAT(hit.first, Catch::Matchers::WithinAbs(2.0, 1e-6)); - origin = {0.0, 0.0, 0.0}; - REQUIRE(xdg->point_in_volume(volume, origin)); + origin = {-10.0, 0.0, 0.0}; + hit = xdg->ray_fire(volume, origin, dir); + REQUIRE_THAT(hit.first, Catch::Matchers::WithinAbs(15.0, 1e-6)); + origin = {0.0, 0.0, 0.0}; + REQUIRE(xdg->point_in_volume(volume, origin)); + } } TEST_CASE("MOAB Element Types") @@ -128,133 +138,144 @@ TEST_CASE("MOAB Element Types") TEST_CASE("MOAB Get Surface Mesh") { - std::shared_ptr xdg = XDG::create(MeshLibrary::MOAB); - REQUIRE(xdg->mesh_manager()->mesh_library() == MeshLibrary::MOAB); - const auto& mesh_manager = xdg->mesh_manager(); - mesh_manager->load_file("overlap-edge.h5m"); - mesh_manager->init(); + // Generate one test run per enabled backend + auto rt_backend = GENERATE(RTLibrary::EMBREE, RTLibrary::GPRT); + + DYNAMIC_SECTION(fmt::format("Backend = {}", rt_backend)) { + check_ray_tracer_supported(rt_backend); // skip if backend not enabled at configuration time + std::shared_ptr xdg = XDG::create(MeshLibrary::MOAB, rt_backend); + REQUIRE(xdg->mesh_manager()->mesh_library() == MeshLibrary::MOAB); + const auto& mesh_manager = xdg->mesh_manager(); + mesh_manager->load_file("overlap-edge.h5m"); + mesh_manager->init(); + + float fpTol = 1e-5; + + // Define the expected connectivity and vertices for each surface + std::vector> expected_connectivity = { + {2, 3, 5, 3, 0, 4, 5, 4, 1, 3, 4, 5}, /* Surface 1 */ + {5, 4, 1, 4, 7, 3, 2, 3, 6, 4, 5, 7, 5, 0, 7, 7, 6, 3}, /* Surface 2 */ + {5, 0, 7, 5, 4, 2, 1, 3, 6, 3, 4, 7, 4, 5, 7, 3, 7, 6}, /* Surface 3 */ + {5, 4, 0, 3, 4, 6, 1, 3, 6, 4, 5, 7, 5, 2, 7, 4, 7, 6}, /* Surface 4 */ + {3, 0, 4, 5, 1, 3, 5, 3, 4, 5, 4, 2}, /* Surface 5 */ + {4, 5, 6, 0, 5, 3, 7, 2, 4, 4, 3, 5, 7, 4, 6, 7, 6, 1}, /* Surface 6 */ + {6, 4, 5, 3, 5, 4, 1, 5, 3, 7, 0, 4, 7, 4, 6, 7, 6, 2}, /* Surface 7 */ + {3, 6, 4, 1, 5, 3, 7, 2, 4, 3, 5, 6, 7, 4, 6, 7, 6, 0} /* Surface 8 */ + }; + + std::vector> expected_vertices = { + // Surface 1 + { + {0.5, -0.866025, -1.5}, {-1, 0.0, -1.5}, {0.5, 0.866025, -1.5}, + {0.5, 0, -1.5}, {-0.25, -0.433013, -1.5}, {-0.25, 0.433013, -1.5} + }, + // Surface 2 + { + {0.5, -0.866025, -1.5}, {-1, 0.0, -1.5}, {0, 0, 1.5}, + {-0.333333, 0.0, 0.5}, {-0.666667, 0.0, -0.5}, + {-0.25, -0.433013, -1.5}, {0.166667, -0.288675, 0.5}, {0.333333, -0.57735, -0.5} + }, + // Surface 3 + { + {-1, 0.0, -1.5}, {0, 0, 1.5}, {0.5, 0.866025, -1.5}, + {0.166667, 0.288675, 0.5}, {0.333333, 0.57735, -0.5}, {-0.25, 0.433013, -1.5}, + {-0.333333, 0.0, 0.5}, {-0.666667, 0.0, -0.5} + }, + // Surface 4 + { + {0.5, -0.866025, -1.5}, {0, 0, 1.5}, {0.5, 0.866025, -1.5}, + {0.166667, -0.288675, 0.5}, {0.333333, -0.57735, -0.5}, {0.5, 0, -1.5}, + {0.166667, 0.288675, 0.5}, {0.333333, 0.57735, -0.5} + }, + // Surface 5 + { + {0.5, 1.71603, 1.5}, {0.5, -0.0160254, 1.5}, {-1, 0.85, 1.5}, + {0.5, 0.85, 1.5}, {-0.25, 1.28301, 1.5}, {-0.25, 0.416987, 1.5} + }, + // Surface 6 + { + {0, 0.85, -1.5}, {0.5, -0.0160254, 1.5}, {-1, 0.85, 1.5}, + {-0.333333, 0.85, -0.5}, {-0.666667, 0.85, 0.5}, {0.166667, 0.561325, -0.5}, + {0.333333, 0.27265, 0.5}, {-0.25, 0.416987, 1.5} + }, + // Surface 7 + { + {0.5, 1.71603, 1.5}, {0, 0.85, -1.5}, {-1, 0.85, 1.5}, + {0.166667, 1.13868, -0.5}, {0.333333, 1.42735, 0.5}, {-0.333333, 0.85, -0.5}, + {-0.666667, 0.85, 0.5}, {-0.25, 1.28301, 1.5} + }, + // Surface 8 + { + {0.5, 1.71603, 1.5}, {0, 0.85, -1.5}, {0.5, -0.0160254, 1.5}, + {0.166667, 0.561325, -0.5}, {0.333333, 0.27265, 0.5}, {0.166667, 1.13868, -0.5}, + {0.333333, 1.42735, 0.5}, {0.5, 0.85, 1.5} + } + }; - float fpTol = 1e-5; - - // Define the expected connectivity and vertices for each surface - std::vector> expected_connectivity = { - {2, 3, 5, 3, 0, 4, 5, 4, 1, 3, 4, 5}, /* Surface 1 */ - {5, 4, 1, 4, 7, 3, 2, 3, 6, 4, 5, 7, 5, 0, 7, 7, 6, 3}, /* Surface 2 */ - {5, 0, 7, 5, 4, 2, 1, 3, 6, 3, 4, 7, 4, 5, 7, 3, 7, 6}, /* Surface 3 */ - {5, 4, 0, 3, 4, 6, 1, 3, 6, 4, 5, 7, 5, 2, 7, 4, 7, 6}, /* Surface 4 */ - {3, 0, 4, 5, 1, 3, 5, 3, 4, 5, 4, 2}, /* Surface 5 */ - {4, 5, 6, 0, 5, 3, 7, 2, 4, 4, 3, 5, 7, 4, 6, 7, 6, 1}, /* Surface 6 */ - {6, 4, 5, 3, 5, 4, 1, 5, 3, 7, 0, 4, 7, 4, 6, 7, 6, 2}, /* Surface 7 */ - {3, 6, 4, 1, 5, 3, 7, 2, 4, 3, 5, 6, 7, 4, 6, 7, 6, 0} /* Surface 8 */ - }; - - std::vector> expected_vertices = { - // Surface 1 - { - {0.5, -0.866025, -1.5}, {-1, 0.0, -1.5}, {0.5, 0.866025, -1.5}, - {0.5, 0, -1.5}, {-0.25, -0.433013, -1.5}, {-0.25, 0.433013, -1.5} - }, - // Surface 2 - { - {0.5, -0.866025, -1.5}, {-1, 0.0, -1.5}, {0, 0, 1.5}, - {-0.333333, 0.0, 0.5}, {-0.666667, 0.0, -0.5}, - {-0.25, -0.433013, -1.5}, {0.166667, -0.288675, 0.5}, {0.333333, -0.57735, -0.5} - }, - // Surface 3 - { - {-1, 0.0, -1.5}, {0, 0, 1.5}, {0.5, 0.866025, -1.5}, - {0.166667, 0.288675, 0.5}, {0.333333, 0.57735, -0.5}, {-0.25, 0.433013, -1.5}, - {-0.333333, 0.0, 0.5}, {-0.666667, 0.0, -0.5} - }, - // Surface 4 - { - {0.5, -0.866025, -1.5}, {0, 0, 1.5}, {0.5, 0.866025, -1.5}, - {0.166667, -0.288675, 0.5}, {0.333333, -0.57735, -0.5}, {0.5, 0, -1.5}, - {0.166667, 0.288675, 0.5}, {0.333333, 0.57735, -0.5} - }, - // Surface 5 - { - {0.5, 1.71603, 1.5}, {0.5, -0.0160254, 1.5}, {-1, 0.85, 1.5}, - {0.5, 0.85, 1.5}, {-0.25, 1.28301, 1.5}, {-0.25, 0.416987, 1.5} - }, - // Surface 6 - { - {0, 0.85, -1.5}, {0.5, -0.0160254, 1.5}, {-1, 0.85, 1.5}, - {-0.333333, 0.85, -0.5}, {-0.666667, 0.85, 0.5}, {0.166667, 0.561325, -0.5}, - {0.333333, 0.27265, 0.5}, {-0.25, 0.416987, 1.5} - }, - // Surface 7 - { - {0.5, 1.71603, 1.5}, {0, 0.85, -1.5}, {-1, 0.85, 1.5}, - {0.166667, 1.13868, -0.5}, {0.333333, 1.42735, 0.5}, {-0.333333, 0.85, -0.5}, - {-0.666667, 0.85, 0.5}, {-0.25, 1.28301, 1.5} - }, - // Surface 8 - { - {0.5, 1.71603, 1.5}, {0, 0.85, -1.5}, {0.5, -0.0160254, 1.5}, - {0.166667, 0.561325, -0.5}, {0.333333, 0.27265, 0.5}, {0.166667, 1.13868, -0.5}, - {0.333333, 1.42735, 0.5}, {0.5, 0.85, 1.5} + size_t surface_index = 0; + for (const auto surface : mesh_manager->surfaces()) { + auto surfaceMesh = mesh_manager->get_surface_mesh(surface); + auto vertices = surfaceMesh.first; + auto connectivity = surfaceMesh.second; + + // Test connectivity + REQUIRE(connectivity.size() == expected_connectivity[surface_index].size()); + for (size_t i = 0; i < connectivity.size(); ++i) { + REQUIRE(connectivity[i] == expected_connectivity[surface_index][i]); } - }; - - size_t surface_index = 0; - for (const auto surface : mesh_manager->surfaces()) { - auto surfaceMesh = mesh_manager->get_surface_mesh(surface); - auto vertices = surfaceMesh.first; - auto connectivity = surfaceMesh.second; - - // Test connectivity - REQUIRE(connectivity.size() == expected_connectivity[surface_index].size()); - for (size_t i = 0; i < connectivity.size(); ++i) { - REQUIRE(connectivity[i] == expected_connectivity[surface_index][i]); - } - // Test vertices - REQUIRE(vertices.size() == expected_vertices[surface_index].size()); - for (size_t i = 0; i < vertices.size(); ++i) { - REQUIRE_THAT(vertices[i].x, Catch::Matchers::WithinAbs(expected_vertices[surface_index][i].x, fpTol)); - REQUIRE_THAT(vertices[i].y, Catch::Matchers::WithinAbs(expected_vertices[surface_index][i].y, fpTol)); - REQUIRE_THAT(vertices[i].z, Catch::Matchers::WithinAbs(expected_vertices[surface_index][i].z, fpTol)); - } + // Test vertices + REQUIRE(vertices.size() == expected_vertices[surface_index].size()); + for (size_t i = 0; i < vertices.size(); ++i) { + REQUIRE_THAT(vertices[i].x, Catch::Matchers::WithinAbs(expected_vertices[surface_index][i].x, fpTol)); + REQUIRE_THAT(vertices[i].y, Catch::Matchers::WithinAbs(expected_vertices[surface_index][i].y, fpTol)); + REQUIRE_THAT(vertices[i].z, Catch::Matchers::WithinAbs(expected_vertices[surface_index][i].z, fpTol)); + } - ++surface_index; + ++surface_index; + } } } TEST_CASE("TEST MOAB Find Element Method") { - std::shared_ptr xdg = XDG::create(MeshLibrary::MOAB); - REQUIRE(xdg->mesh_manager()->mesh_library() == MeshLibrary::MOAB); - const auto& mesh_manager = xdg->mesh_manager(); - mesh_manager->load_file("jezebel.h5m"); - mesh_manager->init(); - xdg->prepare_raytracer(); - - REQUIRE(mesh_manager->num_volume_elements() == 10333); + // Generate one test run per enabled backend + auto rt_backend = GENERATE(RTLibrary::EMBREE); // TODO add GPRT once find element is implemented with GPRT + + DYNAMIC_SECTION(fmt::format("Backend = {}", rt_backend)) { + check_ray_tracer_supported(rt_backend); // skip if backend not enabled at configuration time + std::shared_ptr xdg = XDG::create(MeshLibrary::MOAB, RTLibrary::EMBREE); + REQUIRE(xdg->ray_tracing_interface()->library() == RTLibrary::EMBREE); + REQUIRE(xdg->mesh_manager()->mesh_library() == MeshLibrary::MOAB); + const auto& mesh_manager = xdg->mesh_manager(); + mesh_manager->load_file("jezebel.h5m"); + mesh_manager->init(); + xdg->prepare_raytracer(); MeshID volume = 1; REQUIRE(mesh_manager->num_volume_elements(1) == 10333); - MeshID element = xdg->find_element(volume, {0.0, 0.0, 100.0}); - REQUIRE(element == ID_NONE); // should not find an element since the point is outside the volume - - element = xdg->find_element(volume, {0.0, 0.0, 0.0}); - REQUIRE(element != ID_NONE); // should find an element - - // test the next_element method - auto next_element = xdg->mesh_manager()->next_element(element, {0.0, 0.0, 0.0}, {0.0, 0.0, 1.0}); - REQUIRE(next_element.first != ID_NONE); - REQUIRE(next_element.second != INFTY); - - // test the walk_elements method - auto walk_elements = xdg->mesh_manager()->walk_elements(element, {0.0, 0.0, 0.0}, {0.0, 0.0, 1.0}, 100.0); - // get the sum of the distances - double distance = std::accumulate(walk_elements.begin(), walk_elements.end(), 0.0, - [](double total, const auto& segment) { return total + segment.second; }); - REQUIRE(distance > 0.0); - REQUIRE(distance <= 100.0); - for (const auto& segment : walk_elements) { - REQUIRE(segment.first != ID_NONE); - REQUIRE(segment.second >= 0.0); + MeshID element = xdg->find_element(volume, {0.0, 0.0, 100.0}); + REQUIRE(element == ID_NONE); // should not find an element since the point is outside the volume + + element = xdg->find_element(volume, {0.0, 0.0, 0.0}); + REQUIRE(element != ID_NONE); // should find an element + + // test the next_element method + auto next_element = xdg->mesh_manager()->next_element(element, {0.0, 0.0, 0.0}, {0.0, 0.0, 1.0}); + REQUIRE(next_element.first != ID_NONE); + REQUIRE(next_element.second != INFTY); + + // test the walk_elements method + auto walk_elements = xdg->mesh_manager()->walk_elements(element, {0.0, 0.0, 0.0}, {0.0, 0.0, 1.0}, 100.0); + // get the sum of the distances + double distance = std::accumulate(walk_elements.begin(), walk_elements.end(), 0.0, + [](double total, const auto& segment) { return total + segment.second; }); + REQUIRE(distance > 0.0); + REQUIRE(distance <= 100.0); + for (const auto& segment : walk_elements) { + REQUIRE(segment.first != ID_NONE); + REQUIRE(segment.second >= 0.0); + } } } \ No newline at end of file diff --git a/tests/test_point_in_volume.cpp b/tests/test_point_in_volume.cpp index 318db1e9..a5e75b6a 100644 --- a/tests/test_point_in_volume.cpp +++ b/tests/test_point_in_volume.cpp @@ -1,62 +1,78 @@ // for testing #include +#include + // xdg includes +#include "xdg/constants.h" #include "xdg/mesh_manager_interface.h" -#include "xdg/ray_tracing_interface.h" -#include "xdg/embree/ray_tracer.h" +#include "util.h" #include "mesh_mock.h" using namespace xdg; -TEST_CASE("Test Point in Volume") +// ---------- single test, sections per backend -------------------------------- + +TEST_CASE("Point-in-volume on MeshMock", "[piv][mock]") { - std::shared_ptr mm = std::make_shared(false); - mm->init(); // this should do nothing, just good practice to call it - REQUIRE(mm->mesh_library() == MeshLibrary::MOCK); - - std::shared_ptr rti = std::make_shared(); - auto [volume_tree, element_tree] = rti->register_volume(mm, mm->volumes()[0]); - REQUIRE(volume_tree != ID_NONE); - REQUIRE(element_tree == ID_NONE); - - Position point {0.0, 0.0, 0.0}; - bool result = rti->point_in_volume(volume_tree, point); - REQUIRE(result == true); - - point = {0.0, 0.0, 1000.0}; - result = rti->point_in_volume(volume_tree, point); - REQUIRE(result == false); - - // test a point just inside the positive x boundary - point = {4.0 - 1e-06, 0.0, 0.0}; - result = rti->point_in_volume(volume_tree, point); - REQUIRE(result == true); - - // test a point just outside on the positive x boundary - // no direction - point = {5.001, 0.0, 0,0}; - result = rti->point_in_volume(volume_tree, point); - REQUIRE(result == false); - - // test a point on the positive x boundary - // and provide a direction - point = {5.0, 0.0, 0.0}; - Direction dir = {1.0, 0.0, 0.0}; - result = rti->point_in_volume(volume_tree, point, &dir); - REQUIRE(result == true); - - // test a point just outside the positive x boundary - // and provide a direction - point = {5.1, 0.0, 0.0}; - dir = {1.0, 0.0, 0.0}; - result = rti->point_in_volume(volume_tree, point, &dir); - REQUIRE(result == false); - - // test a point just outside the positive x boundary, - // flip the direction - point = {5.1, 0.0, 0.0}; - dir = {-1.0, 0.0, 0.0}; - result = rti->point_in_volume(volume_tree, point, &dir); - REQUIRE(result == false); + // Generate one test run per enabled backend + auto rt_backend = GENERATE(RTLibrary::EMBREE, RTLibrary::GPRT); + + DYNAMIC_SECTION(fmt::format("Backend = {}", rt_backend)) { + check_ray_tracer_supported(rt_backend); // skip if backend not enabled at configuration time + auto rti = create_raytracer(rt_backend); + REQUIRE(rti); + rti->init(); + + // Keep MeshMock usage consistent across backends + auto mm = std::make_shared(false); + mm->init(); + REQUIRE(mm->mesh_library() == MeshLibrary::MOCK); + + auto [volume_tree, element_tree] = rti->register_volume(mm, mm->volumes()[0]); + REQUIRE(volume_tree != ID_NONE); + REQUIRE(element_tree == ID_NONE); + + rti->init(); // Ensure ray tracer is initialized (e.g. build SBT for GPRT) + + Position point {0.0, 0.0, 0.0}; + bool result = rti->point_in_volume(volume_tree, point); + REQUIRE(result == true); + + point = {0.0, 0.0, 1000.0}; + result = rti->point_in_volume(volume_tree, point); + REQUIRE(result == false); + + // test a point just inside the positive x boundary + point = {4.0 - 1e-6, 0.0, 0.0}; + result = rti->point_in_volume(volume_tree, point); + REQUIRE(result == true); + + // test a point just outside on the positive x boundary + // no direction + point = {5.001, 0.0, 0.0}; + result = rti->point_in_volume(volume_tree, point); + REQUIRE(result == false); + + // test a point on the positive x boundary + // and provide a direction + point = {5.0, 0.0, 0.0}; + Direction dir {1.0, 0.0, 0.0}; + result = rti->point_in_volume(volume_tree, point, &dir); + REQUIRE(result == true); + + // test a point just outside the positive x boundary + // and provide a direction + point = {5.1, 0.0, 0.0}; + dir = {1.0, 0.0, 0.0}; + result = rti->point_in_volume(volume_tree, point, &dir); + REQUIRE(result == false); + + // test a point just outside the positive x boundary, + // flip the direction + point = {5.1, 0.0, 0.0}; + dir = {-1.0, 0.0, 0.0}; + result = rti->point_in_volume(volume_tree, point, &dir); + REQUIRE(result == false); + } } \ No newline at end of file diff --git a/tests/test_ray_fire.cpp b/tests/test_ray_fire.cpp index 5fd57f02..8c360b87 100644 --- a/tests/test_ray_fire.cpp +++ b/tests/test_ray_fire.cpp @@ -1,102 +1,106 @@ - // for testing #include #include +#include + // xdg includes #include "xdg/constants.h" #include "xdg/mesh_manager_interface.h" -#include "xdg/embree/ray_tracer.h" - #include "mesh_mock.h" +#include "util.h" using namespace xdg; -TEST_CASE("Test Ray Fire Mesh Mock") -{ - std::shared_ptr mm = std::make_shared(false); - mm->init(); // this should do nothing, just good practice to call it - REQUIRE(mm->mesh_library() == MeshLibrary::MOCK); - - std::shared_ptr rti = std::make_shared(); - auto [volume_tree, element_tree] = rti->register_volume(mm, mm->volumes()[0]); - REQUIRE(volume_tree != ID_NONE); - REQUIRE(element_tree == ID_NONE); - - Position origin {0.0, 0.0, 0.0}; - Direction direction {1.0, 0.0, 0.0}; - std::pair intersection; - - // fire from the origin toward each face, ensuring that the intersection distances are correct - intersection = rti->ray_fire(volume_tree, origin, direction); - REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(5.0, 1e-6)); - - direction *= -1; - intersection = rti->ray_fire(volume_tree, origin, direction); - REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(2.0, 1e-6)); - - direction = {0.0, 1.0, 0.0}; - intersection = rti->ray_fire(volume_tree, origin, direction); - REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(6.0, 1e-6)); - - direction *= -1; - intersection = rti->ray_fire(volume_tree, origin, direction); - REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(3.0, 1e-6)); - - direction = {0.0, 0.0, 1.0}; - intersection = rti->ray_fire(volume_tree, origin, direction); - REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(7.0, 1e-6)); - - direction *= -1; - intersection = rti->ray_fire(volume_tree, origin, direction); - REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(4.0, 1e-6)); - - // fire from the outside of the cube toward each face, ensuring that the intersection distances are correct - // rays should skip entering intersections and intersect with the far side of the cube - origin = {-10.0, 0.0, 0.0}; - direction = {1.0, 0.0, 0.0}; - intersection = rti->ray_fire(volume_tree, origin, direction); - REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(15.0, 1e-6)); - - origin = {10.0, 0.0, 0.0}; - direction = {-1.0, 0.0, 0.0}; - intersection = rti->ray_fire(volume_tree, origin, direction); - REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(12.0, 1e-6)); - - // fire from the outside of the cube toward each face, ensuring that the intersection distances are correct - // in this case rays are fired with a HitOrientation::ENTERING. Rays should hit the first surface intersected - origin = {-10.0, 0.0, 0.0}; - direction = {1.0, 0.0, 0.0}; - intersection = rti->ray_fire(volume_tree, origin, direction, INFTY, HitOrientation::ENTERING); - REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(8.0, 1e-6)); - - origin = {10.0, 0.0, 0.0}; - direction = {-1.0, 0.0, 0.0}; - intersection = rti->ray_fire(volume_tree, origin, direction, INFTY, HitOrientation::ENTERING); - REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(5.0, 1e-6)); - - // limit distance of the ray, shouldn't get a hit - origin = {0.0, 0.0, 0.0}; - direction = {1.0, 0.0, 0.0}; - intersection = rti->ray_fire(volume_tree, origin, direction, 4.5); - REQUIRE(intersection.second == ID_NONE); - - // if the distance is just enough, we should still get a hit - // limit distance of the ray, shouldn't get a hit - origin = {0.0, 0.0, 0.0}; - direction = {1.0, 0.0, 0.0}; - intersection = rti->ray_fire(volume_tree, origin, direction, 5.1); - REQUIRE(intersection.second != ID_NONE); - REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(5.0, 1e-6)); - - // Test excluding primitives, fire a ray from the origin and log the hit face - // By providing the hit face as an excluded primitive in a subsequent ray fire, - // there should be no intersection returned - std::vector exclude_primitives; - intersection = rti->ray_fire(volume_tree, origin, direction, INFTY, HitOrientation::EXITING, &exclude_primitives); - REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(5.0, 1e-6)); - REQUIRE(exclude_primitives.size() == 1); - - intersection = rti->ray_fire(volume_tree, origin, direction, INFTY, HitOrientation::EXITING, &exclude_primitives); - REQUIRE(intersection.second == ID_NONE); -} +// ------- single test, multiple sections (one per built backend) -------------- + +TEST_CASE("Ray Fire on MeshMock (per-backend sections)", "[rayfire][mock]") { + // Generate one test run per enabled backend + auto rt_backend = GENERATE(RTLibrary::EMBREE, RTLibrary::GPRT); + + DYNAMIC_SECTION(fmt::format("Backend = {}", rt_backend)) { + check_ray_tracer_supported(rt_backend); // skip if backend not enabled at configuration time + auto rti = create_raytracer(rt_backend); + REQUIRE(rti); + + auto mm = std::make_shared(false); + mm->init(); + REQUIRE(mm->mesh_library() == MeshLibrary::MOCK); + + auto [volume_tree, element_tree] = rti->register_volume(mm, mm->volumes()[0]); + REQUIRE(volume_tree != ID_NONE); + REQUIRE(element_tree == ID_NONE); + + rti->init(); // Ensure ray tracer is initialized (e.g. build SBT for GPRT) + Position origin {0.0, 0.0, 0.0}; + Direction direction {1.0, 0.0, 0.0}; + std::pair intersection; + + intersection = rti->ray_fire(volume_tree, origin, direction); + REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(5.0, 1e-6)); + + direction *= -1; + intersection = rti->ray_fire(volume_tree, origin, direction); + REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(2.0, 1e-6)); + + direction = {0.0, 1.0, 0.0}; + intersection = rti->ray_fire(volume_tree, origin, direction); + REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(6.0, 1e-6)); + direction *= -1; + intersection = rti->ray_fire(volume_tree, origin, direction); + REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(3.0, 1e-6)); + + direction = {0.0, 0.0, 1.0}; + intersection = rti->ray_fire(volume_tree, origin, direction); + REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(7.0, 1e-6)); + direction *= -1; + intersection = rti->ray_fire(volume_tree, origin, direction); + REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(4.0, 1e-6)); + + // fire from the outside of the cube toward each face, ensuring that the intersection distances are correct + // rays should skip entering intersections and intersect with the far side of the cube + origin = {-10.0, 0.0, 0.0}; + direction = {1.0, 0.0, 0.0}; + intersection = rti->ray_fire(volume_tree, origin, direction); + REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(15.0, 1e-6)); + + origin = {10.0, 0.0, 0.0}; + direction = {-1.0, 0.0, 0.0}; + intersection = rti->ray_fire(volume_tree, origin, direction); + REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(12.0, 1e-6)); + + // fire from the outside of the cube toward each face, ensuring that the intersection distances are correct + // in this case rays are fired with a HitOrientation::ENTERING. Rays should hit the first surface intersected + origin = {-10.0, 0.0, 0.0}; + direction = {1.0, 0.0, 0.0}; + intersection = rti->ray_fire(volume_tree, origin, direction, INFTY, HitOrientation::ENTERING); + REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(8.0, 1e-6)); + + origin = {10.0, 0.0, 0.0}; + direction = {-1.0, 0.0, 0.0}; + intersection = rti->ray_fire(volume_tree, origin, direction, INFTY, HitOrientation::ENTERING); + REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(5.0, 1e-6)); + + // limit distance of the ray, shouldn't get a hit + origin = {0.0, 0.0, 0.0}; + direction = {1.0, 0.0, 0.0}; + intersection = rti->ray_fire(volume_tree, origin, direction, 4.5); + REQUIRE(intersection.second == ID_NONE); + + // if the distance is just enough, we should still get a hit + intersection = rti->ray_fire(volume_tree, origin, direction, 5.1); + REQUIRE(intersection.second != ID_NONE); + REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(5.0, 1e-6)); + + // Test excluding primitives, fire a ray from the origin and log the hit face + // By providing the hit face as an excluded primitive in a subsequent ray fire, + // there should be no intersection returned + std::vector exclude_primitives; + intersection = rti->ray_fire(volume_tree, origin, direction, INFTY, HitOrientation::EXITING, &exclude_primitives); + REQUIRE_THAT(intersection.first, Catch::Matchers::WithinAbs(5.0, 1e-6)); + REQUIRE(exclude_primitives.size() == 1); + + intersection = rti->ray_fire(volume_tree, origin, direction, INFTY, HitOrientation::EXITING, &exclude_primitives); + REQUIRE(intersection.second == ID_NONE); + } +} \ No newline at end of file diff --git a/tests/util.h b/tests/util.h index 1256212d..2d5f6bc0 100644 --- a/tests/util.h +++ b/tests/util.h @@ -1,5 +1,11 @@ #include +#include + +#include "xdg/constants.h" +#include "xdg/ray_tracers.h" +#include "vulkan_probe.h" + static std::random_device rd; static std::mt19937 gen(rd()); @@ -9,3 +15,33 @@ inline double rand_double(double min, double max) return dis(gen); } +inline void check_ray_tracer_supported(xdg::RTLibrary rt) { + #ifndef XDG_ENABLE_EMBREE + if (rt == xdg::RTLibrary::EMBREE) { + SKIP("XDG not built with Embree backend; skipping Embree tests."); + } + #endif + + #ifndef XDG_ENABLE_GPRT + if (rt == xdg::RTLibrary::GPRT) { + SKIP("XDG not built with GPRT backend; skipping GPRT tests."); + } + #else // XDG_ENABLE_GPRT + if (rt == xdg::RTLibrary::GPRT && !system_has_vk_device()) { + SKIP("No Vulkan device found; skipping GPRT tests."); + } + #endif +} + +// Factory function to create ray tracer based on which library selected +inline std::shared_ptr create_raytracer(xdg::RTLibrary rt) { + #ifdef XDG_ENABLE_EMBREE + if (rt == xdg::RTLibrary::EMBREE) + return std::make_shared(); + #endif + + #ifdef XDG_ENABLE_GPRT + if (rt == xdg::RTLibrary::GPRT) + return std::make_shared(); + #endif +} \ No newline at end of file diff --git a/tests/vulkan_probe.h b/tests/vulkan_probe.h new file mode 100644 index 00000000..cfc6ff44 --- /dev/null +++ b/tests/vulkan_probe.h @@ -0,0 +1,41 @@ +#pragma once + +#ifdef XDG_ENABLE_GPRT +#include +#include + +inline bool system_has_vk_device(uint32_t min_instance = VK_API_VERSION_1_1) { + uint32_t loaderVer = VK_API_VERSION_1_0; + vkEnumerateInstanceVersion(&loaderVer); + + // Create VK instance + VkApplicationInfo app{ VK_STRUCTURE_TYPE_APPLICATION_INFO }; + app.pApplicationName = "vk-probe"; + app.applicationVersion = 1; + app.pEngineName = "probe"; + app.engineVersion = 1; + app.apiVersion = VK_API_VERSION_1_0; + + VkInstanceCreateInfo ici{ VK_STRUCTURE_TYPE_INSTANCE_CREATE_INFO }; + ici.pApplicationInfo = &app; + + VkInstance instance = VK_NULL_HANDLE; + if (vkCreateInstance(&ici, nullptr, &instance) != VK_SUCCESS || !instance) + return false; + + // Look for at least one physical device + uint32_t devCount = 0; + VkResult r = vkEnumeratePhysicalDevices(instance, &devCount, nullptr); + + // Clean up the instance before returning + vkDestroyInstance(instance, nullptr); + + // Check for errors and non-zero device count + if (r != VK_SUCCESS || devCount == 0) { + return false; + } + + return true; // VK device detected +} + +#endif diff --git a/tools/CMakeLists.txt b/tools/CMakeLists.txt index 514f536f..121e2521 100644 --- a/tools/CMakeLists.txt +++ b/tools/CMakeLists.txt @@ -1,6 +1,3 @@ - - - set(TOOL_NAMES particle_sim ray_fire diff --git a/tools/particle_sim.cpp b/tools/particle_sim.cpp index 22c4be8d..e5e7dc7e 100644 --- a/tools/particle_sim.cpp +++ b/tools/particle_sim.cpp @@ -30,18 +30,21 @@ args.add_argument("-m", "--mfp") .default_value(1.0) .help("Mean free path of the particles").scan<'g', double>(); -args.add_argument("-l", "--library") +args.add_argument("-m", "--mesh-library") .help("Mesh library to use. One of (MOAB, LIBMESH)") .default_value("MOAB"); - try { - args.parse_args(argc, argv); - } - catch (const std::runtime_error& err) { - std::cout << err.what() << std::endl; - std::cout << args; - exit(0); - } +args.add_argument("-r", "--rt-library") + .help("Ray tracing library to use. One of (EMBREE, GPRT)") + .default_value("EMBREE"); +try { + args.parse_args(argc, argv); +} +catch (const std::runtime_error& err) { + std::cout << err.what() << std::endl; + std::cout << args; + exit(0); +} // Problem Setup srand48(42); @@ -49,22 +52,38 @@ srand48(42); SimulationData sim_data; // create a mesh manager -std::shared_ptr xdg {nullptr}; -if (args.get("--library") == "MOAB") - xdg = XDG::create(MeshLibrary::MOAB); -else if (args.get("--library") == "LIBMESH") - xdg = XDG::create(MeshLibrary::LIBMESH); +std::string mesh_str = args.get("--mesh-library"); +std::string rt_str = args.get("--rt-library"); + +RTLibrary rt_lib; +if (rt_str == "EMBREE") + rt_lib = RTLibrary::EMBREE; +else if (rt_str == "GPRT") + rt_lib = RTLibrary::GPRT; else - fatal_error("Invalid mesh library {} specified", args.get("--library")); - -sim_data.xdg_ = xdg; + fatal_error("Invalid ray tracing library '{}' specified", rt_str); + +MeshLibrary mesh_lib; +if (mesh_str == "MOAB") + mesh_lib = MeshLibrary::MOAB; +else if (mesh_str == "LIBMESH") { + mesh_lib = MeshLibrary::LIBMESH; + if (rt_lib == RTLibrary::GPRT) + fatal_error("LibMesh is not currently supported with GPRT"); +} +else + fatal_error("Invalid mesh library '{}' specified", mesh_str); +// create an XDG instance with the specified mesh and ray tracing library +std::shared_ptr xdg = XDG::create(mesh_lib, rt_lib); const auto& mm = xdg->mesh_manager(); mm->load_file(args.get("filename")); mm->init(); mm->parse_metadata(); xdg->prepare_raytracer(); +sim_data.xdg_ = xdg; + // update the mean free path sim_data.mfp_ = args.get("--mfp"); diff --git a/tools/point_in_volume.cpp b/tools/point_in_volume.cpp index abbdddd2..e32ee9f7 100644 --- a/tools/point_in_volume.cpp +++ b/tools/point_in_volume.cpp @@ -14,7 +14,6 @@ using namespace xdg; int main(int argc, char** argv) { - argparse::ArgumentParser args("XDG Point in Volume Tool", "1.0", argparse::default_arguments::help); args.add_argument("filename") @@ -28,14 +27,22 @@ int main(int argc, char** argv) { .implicit_value(true) .help("List all volumes in the file and exit"); - args.add_argument("-p", "--position") + args.add_argument("-o", "-p", "--origin", "--position") .default_value(std::vector{0.0, 0.0, 0.0}) - .help("Ray origin").scan<'g', double>().nargs(3); + .help("Ray origin/position").scan<'g', double>().nargs(3); args.add_argument("-d", "--direction") .default_value(std::vector{0.0, 0.0, 1.0}) .help("Ray direction").scan<'g', double>().nargs(3); + args.add_argument("-m", "--mesh-library") + .help("Mesh library to use. One of (MOAB, LIBMESH)") + .default_value("MOAB"); + + args.add_argument("-r", "--rt-library") + .help("Ray tracing library to use. One of (EMBREE, GPRT)") + .default_value("EMBREE"); + try { args.parse_args(argc, argv); } @@ -44,14 +51,37 @@ int main(int argc, char** argv) { std::cout << args; exit(0); } + +std::string mesh_str = args.get("--mesh-library"); +std::string rt_str = args.get("--rt-library"); + +RTLibrary rt_lib; +if (rt_str == "EMBREE") + rt_lib = RTLibrary::EMBREE; +else if (rt_str == "GPRT") + rt_lib = RTLibrary::GPRT; +else + fatal_error("Invalid ray tracing library '{}' specified", rt_str); + +MeshLibrary mesh_lib; +if (mesh_str == "MOAB") + mesh_lib = MeshLibrary::MOAB; +else if (mesh_str == "LIBMESH") { + mesh_lib = MeshLibrary::LIBMESH; + if (rt_lib == RTLibrary::GPRT) + fatal_error("LibMesh is not currently supported with GPRT"); +} +else + fatal_error("Invalid mesh library '{}' specified", mesh_str); // create a mesh manager - std::shared_ptr xdg = XDG::create(MeshLibrary::MOAB); + std::shared_ptr xdg = XDG::create(mesh_lib, rt_lib); const auto& mm = xdg->mesh_manager(); mm->load_file(args.get("filename")); mm->init(); mm->parse_metadata(); - xdg->prepare_raytracer(); + + auto rti = xdg->ray_tracing_interface(); if (args.get("--list")) { std::cout << "Volumes: " << std::endl; @@ -65,6 +95,8 @@ int main(int argc, char** argv) { Position position = args.get>("--position"); Direction direction = args.get>("--direction"); + xdg->prepare_volume_for_raytracing(volume); + if (xdg->point_in_volume(volume, position, &direction)) { std::cout << "Point " << position << " is in Volume " << volume << " (True)" << std::endl; } else { diff --git a/tools/ray_fire.cpp b/tools/ray_fire.cpp index 6f9879a2..64ca8f9f 100644 --- a/tools/ray_fire.cpp +++ b/tools/ray_fire.cpp @@ -2,6 +2,7 @@ #include #include #include +#include #include "xdg/error.h" #include "xdg/mesh_manager_interface.h" @@ -36,6 +37,14 @@ int main(int argc, char** argv) { .default_value(std::vector{0.0, 0.0, 1.0}) .help("Ray direction").scan<'g', double>().nargs(3); + args.add_argument("-m", "--mesh-library") + .help("Mesh library to use. One of (MOAB, LIBMESH)") + .default_value("MOAB"); + + args.add_argument("-r", "--rt-library") + .help("Ray tracing library to use. One of (EMBREE, GPRT)") + .default_value("EMBREE"); + try { args.parse_args(argc, argv); } @@ -44,15 +53,38 @@ int main(int argc, char** argv) { std::cout << args; exit(0); } - + +std::string mesh_str = args.get("--mesh-library"); +std::string rt_str = args.get("--rt-library"); + +RTLibrary rt_lib; +if (rt_str == "EMBREE") + rt_lib = RTLibrary::EMBREE; +else if (rt_str == "GPRT") + rt_lib = RTLibrary::GPRT; +else + fatal_error("Invalid ray tracing library '{}' specified", rt_str); + +MeshLibrary mesh_lib; +if (mesh_str == "MOAB") + mesh_lib = MeshLibrary::MOAB; +else if (mesh_str == "LIBMESH") { + mesh_lib = MeshLibrary::LIBMESH; + if (rt_lib == RTLibrary::GPRT) + fatal_error("LibMesh is not currently supported with GPRT"); +} +else + fatal_error("Invalid mesh library '{}' specified", mesh_str); // create a mesh manager - std::shared_ptr xdg = XDG::create(MeshLibrary::MOAB); + std::shared_ptr xdg = XDG::create(mesh_lib, rt_lib); const auto& mm = xdg->mesh_manager(); mm->load_file(args.get("filename")); mm->init(); mm->parse_metadata(); + auto rti = xdg->ray_tracing_interface(); + if (args.get("--list")) { std::cout << "Volumes: " << std::endl; for (auto volume : mm->volumes()) { @@ -73,7 +105,7 @@ int main(int argc, char** argv) { auto result = xdg->ray_fire(volume, origin, direction); - std::cout << "Distance: " << result.first << std::endl; + std::cout << std::setprecision(17) << "Distance: " << result.first << std::endl; std::cout << "Surface: " << result.second << std::endl; return 0; diff --git a/vendor/GPRT b/vendor/GPRT new file mode 160000 index 00000000..f1e95e41 --- /dev/null +++ b/vendor/GPRT @@ -0,0 +1 @@ +Subproject commit f1e95e4188cde591547d6b4a33a70bf2afaeec59