From b6372072e76b98ae9f0c128a3c372f959635d5dd Mon Sep 17 00:00:00 2001 From: reiter Date: Tue, 3 Feb 2026 13:11:53 +0100 Subject: [PATCH 01/10] Use 3D direction in 2D angle calculation --- examples/disk2D/disk2D.cpp | 1 + include/viennaray/rayBoundary.hpp | 28 ++++++----- include/viennaray/raySourceRandom.hpp | 26 ++++------ include/viennaray/rayTraceDisk.hpp | 17 ++++--- include/viennaray/rayTraceKernel.hpp | 55 ++++++++++----------- include/viennaray/rayUtil.hpp | 69 ++++++++++++++------------- 6 files changed, 95 insertions(+), 101 deletions(-) diff --git a/examples/disk2D/disk2D.cpp b/examples/disk2D/disk2D.cpp index 48517f7..8fda7c6 100644 --- a/examples/disk2D/disk2D.cpp +++ b/examples/disk2D/disk2D.cpp @@ -60,6 +60,7 @@ int main() { // Extract the normalized hit counts for each geometry point auto &flux = rayTracer.getLocalData().getVectorData("flux"); rayTracer.normalizeFlux(flux, NormalizationType::SOURCE); + rayTracer.smoothFlux(flux, 1); rayInternal::writeVTK("trenchResult.vtk", points, flux); diff --git a/include/viennaray/rayBoundary.hpp b/include/viennaray/rayBoundary.hpp index b0e5605..83e2927 100644 --- a/include/viennaray/rayBoundary.hpp +++ b/include/viennaray/rayBoundary.hpp @@ -70,6 +70,7 @@ template class Boundary { assert(false && "Correctness Assumption"); } else { + // D == 3 if (primID <= 3) { if (boundaryConds_[0] == BoundaryCondition::REFLECTIVE_BOUNDARY) { // use specular reflection @@ -151,11 +152,11 @@ template class Boundary { private: static Vec3D getNewOrigin(const RTCRay &ray) { - assert(IsNormalized(Vec3D{ray.dir_x, ray.dir_y, ray.dir_z}) && - "MetaGeometry: direction not normalized"); - auto xx = ray.org_x + ray.dir_x * ray.tfar; - auto yy = ray.org_y + ray.dir_y * ray.tfar; - auto zz = ray.org_z + ray.dir_z * ray.tfar; + auto dir = *reinterpret_cast *>(&ray.dir_x); + assert(IsNormalized(dir) && "Ray direction not normalized"); + auto xx = ray.org_x + dir[0] * ray.tfar; + auto yy = ray.org_y + dir[1] * ray.tfar; + auto zz = ray.org_z + dir[2] * ray.tfar; return Vec3D{xx, yy, zz}; } @@ -257,16 +258,17 @@ template class Boundary { } static void reflectRay(RTCRayHit &rayHit) { - auto dir = *reinterpret_cast *>( - &rayHit.ray.dir_x); - auto normal = *reinterpret_cast *>( - &rayHit.hit.Ng_x); - Normalize(dir); + auto dir = *reinterpret_cast *>(&rayHit.ray.dir_x); + assert(IsNormalized(dir) && "Direction not normalized"); + assert(D == 3 || dir[2] == 0.f && "In 2D case, z-component must be 0"); + auto normal = *reinterpret_cast *>(&rayHit.hit.Ng_x); + // Normalize(dir); Normalize(normal); - auto const direction = - ReflectionSpecular(dir, normal); + auto const direction = ReflectionSpecular(dir, normal); auto const origin = getNewOrigin(rayHit.ray); - rayInternal::fillRayDirection(rayHit.ray, direction); + // use here because z-component should already be 0 in D=2 case + assert(D == 3 || dir[2] == 0.f && "In 2D case, z-component must be 0"); + rayInternal::fillRayDirection<3>(rayHit.ray, direction); rayInternal::fillRayPosition(rayHit.ray, origin); } diff --git a/include/viennaray/raySourceRandom.hpp b/include/viennaray/raySourceRandom.hpp index 2a56e8d..6e2ed8f 100644 --- a/include/viennaray/raySourceRandom.hpp +++ b/include/viennaray/raySourceRandom.hpp @@ -73,17 +73,14 @@ class SourceRandom : public Source { auto r1 = uniDist(rngState); auto r2 = uniDist(rngState); + NumericType sinPhi, cosPhi; + rayInternal::sincos(M_PI * 2. * r1, sinPhi, cosPhi); const NumericType cosTheta = std::pow(r2, ee_); const NumericType sinTheta = std::sqrt(1. - cosTheta * cosTheta); - direction[rayDir_] = posNeg_ * cosTheta; - direction[firstDir_] = std::cos(M_PI * 2. * r1) * sinTheta; - if constexpr (D == 2) { - direction[secondDir_] = 0.; - Normalize(direction); - } else { - direction[secondDir_] = std::sin(M_PI * 2. * r1) * sinTheta; - } + direction[rayDir_] = posNeg_ * cosTheta; + direction[firstDir_] = cosPhi * sinTheta; + direction[secondDir_] = sinPhi * sinTheta; return direction; } @@ -93,15 +90,15 @@ class SourceRandom : public Source { std::uniform_real_distribution uniDist; do { - Vec3D rndDirection{0., 0., 0.}; auto r1 = uniDist(rngState); auto r2 = uniDist(rngState); + NumericType sinPhi, cosPhi; + rayInternal::sincos(M_PI * 2. * r1, sinPhi, cosPhi); const NumericType cosTheta = std::pow(r2, ee_); const NumericType sinTheta = std::sqrt(1. - cosTheta * cosTheta); - rndDirection[0] = cosTheta; - rndDirection[1] = std::cos(M_PI * 2. * r1) * sinTheta; - rndDirection[2] = std::sin(M_PI * 2. * r1) * sinTheta; + Vec3D rndDirection{cosTheta, cosPhi * sinTheta, + sinPhi * sinTheta}; direction[0] = orthonormalBasis_[0][0] * rndDirection[0] + orthonormalBasis_[1][0] * rndDirection[1] + @@ -115,11 +112,6 @@ class SourceRandom : public Source { } while ((posNeg_ < 0. && direction[rayDir_] > 0.) || (posNeg_ > 0. && direction[rayDir_] < 0.)); - if constexpr (D == 2) { - direction[secondDir_] = 0; - Normalize(direction); - } - return direction; } diff --git a/include/viennaray/rayTraceDisk.hpp b/include/viennaray/rayTraceDisk.hpp index cb273c7..38766c0 100644 --- a/include/viennaray/rayTraceDisk.hpp +++ b/include/viennaray/rayTraceDisk.hpp @@ -106,10 +106,9 @@ class TraceDisk final : public Trace { assert(flux.size() == geometry_.getNumPrimitives() && "Unequal number of points in normalizeFlux"); - const auto totalDiskArea = diskRadius_ * diskRadius_ * M_PI; - switch (norm) { case NormalizationType::MAX: { + const auto totalDiskArea = diskRadius_ * diskRadius_ * M_PI; auto maxv = *std::max_element(flux.begin(), flux.end()); #pragma omp parallel for for (int idx = 0; idx < flux.size(); ++idx) { @@ -124,12 +123,12 @@ class TraceDisk final : public Trace { "No source was specified in rayTrace for the normalization."); break; } - NumericType sourceArea = this->pSource_->getSourceArea(); - auto numTotalRays = + const NumericType sourceArea = this->pSource_->getSourceArea(); + const auto numTotalRays = this->config_.numRaysFixed == 0 ? this->pSource_->getNumPoints() * this->config_.numRaysPerPoint : this->config_.numRaysFixed; - NumericType normFactor = sourceArea / numTotalRays; + const NumericType normFactor = sourceArea / numTotalRays; #pragma omp parallel for for (int idx = 0; idx < flux.size(); ++idx) { flux[idx] *= normFactor / geometry_.getDiskArea(idx); @@ -148,6 +147,12 @@ class TraceDisk final : public Trace { int numNeighbors = 1) override { assert(flux.size() == geometry_.getNumPrimitives() && "Unequal number of points in smoothFlux"); + if (numNeighbors < 1) { + VIENNACORE_LOG_DEBUG( + "Number of neighbors for flux smoothing less than 1. Skipping."); + return; + } + auto oldFlux = flux; PointNeighborhood pointNeighborhood; if (numNeighbors == 1) { @@ -171,8 +176,8 @@ class TraceDisk final : public Trace { NumericType vv = oldFlux[idx]; auto const &neighborhood = pointNeighborhood.getNeighborIndices(idx); - NumericType sum = 1.; auto const normal = geometry_.getPrimNormal(idx); + NumericType sum = 1.; for (auto const &nbi : neighborhood) { auto nnormal = geometry_.getPrimNormal(nbi); diff --git a/include/viennaray/rayTraceKernel.hpp b/include/viennaray/rayTraceKernel.hpp index f47700a..7d1f8f7 100644 --- a/include/viennaray/rayTraceKernel.hpp +++ b/include/viennaray/rayTraceKernel.hpp @@ -30,9 +30,8 @@ template class TraceKernel { traceInfo_(traceInfo), dataLog_(dataLog) {} void apply() { + // Create Embree scene auto rtcScene = rtcNewScene(device_); - - // RTC scene flags rtcSetSceneFlags(rtcScene, RTC_SCENE_FLAG_NONE); // Selecting higher build quality results in better rendering performance @@ -49,8 +48,8 @@ template class TraceKernel { size_t geoHits = 0; size_t nonGeoHits = 0; - size_t totalTraces = 0; size_t particleHits = 0; + size_t totalTraces = 0; size_t totalBoundaryHits = 0; size_t totalReflections = 0; size_t raysTerminated = 0; @@ -86,7 +85,7 @@ template class TraceKernel { timer.start(); #pragma omp parallel reduction( \ - + : geoHits, nonGeoHits, totalTraces, particleHits, totalBoundaryHits, \ + + : geoHits, nonGeoHits, particleHits, totalTraces, totalBoundaryHits, \ totalReflections, raysTerminated) shared(threadLocalData) { rtcJoinCommitScene(rtcScene); @@ -124,22 +123,24 @@ template class TraceKernel { // probabilistic weight const NumericType initialRayWeight = pSource_->getInitialRayWeight(idx); NumericType rayWeight = initialRayWeight; + Vec3D rayDirection; unsigned int numReflections = 0; unsigned int boundaryHits = 0; { particle->initNew(rngState); particle->logData(myDataLog); - auto direction = particle->initNewWithDirection(rngState); + rayDirection = particle->initNewWithDirection(rngState); auto originAndDirection = pSource_->getOriginAndDirection(idx, rngState); fillRayPosition(rayHit.ray, originAndDirection[0]); - if (isZero(direction)) { - fillRayDirection(rayHit.ray, originAndDirection[1]); + if (isZero(rayDirection)) { + assert(IsNormalized(originAndDirection[1])); + fillRayDirection(rayHit.ray, originAndDirection[1]); } else { - assert(IsNormalized(direction)); - fillRayDirection(rayHit.ray, direction); + assert(IsNormalized(rayDirection)); + fillRayDirection(rayHit.ray, rayDirection); } } @@ -180,27 +181,22 @@ template class TraceKernel { if (lambda > 0.) { std::uniform_real_distribution dist(0., 1.); NumericType scatterProbability = - 1 - std::exp(-rayHit.ray.tfar / lambda); + 1. - std::exp(-rayHit.ray.tfar / lambda); auto rnd = dist(rngState); if (rnd < scatterProbability) { const auto &ray = rayHit.ray; - const auto origin = Vec3D{ - static_cast(ray.org_x + ray.dir_x * rnd), - static_cast(ray.org_y + ray.dir_y * rnd), - static_cast(ray.org_z + ray.dir_z * rnd)}; - - auto direction = - rayInternal::pickRandomPointOnUnitSphere( + const auto origin = + Vec3Df{static_cast(ray.org_x + ray.dir_x * rnd), + static_cast(ray.org_y + ray.dir_y * rnd), + static_cast(ray.org_z + ray.dir_z * rnd)}; + rayDirection = + rayInternal::pickRandomPointOnUnitSphere( rngState); - if constexpr (D == 2) { - direction[2] = 0.f; - Normalize(direction); - } // Update ray direction and origin fillRayPosition(rayHit.ray, origin); - fillRayDirection(rayHit.ray, direction); + fillRayDirection(rayHit.ray, rayDirection); particleHits++; reflect = true; @@ -224,14 +220,11 @@ template class TraceKernel { const auto hitPoint = Vec3Df{ray.org_x + ray.dir_x * ray.tfar, ray.org_y + ray.dir_y * ray.tfar, ray.org_z + ray.dir_z * ray.tfar}; - - const auto rayDir = - Vec3D{ray.dir_x, ray.dir_y, ray.dir_z}; const auto geomNormal = geometry_.getPrimNormal(rayHit.hit.primID); // Check for backface hit in case of disks if constexpr (geoType == GeometryType::DISK) { - if (DotProduct(rayDir, geomNormal) > 0) { + if (DotProduct(rayDirection, geomNormal) > 0) { // If the dot product of the ray direction and the surface normal // is greater than zero, then we hit the back face of the disk. if (hitFromBack) { @@ -296,21 +289,21 @@ template class TraceKernel { #else auto distRayWeight = rayWeight; #endif - particle->surfaceCollision(distRayWeight, rayDir, normal, + particle->surfaceCollision(distRayWeight, rayDirection, normal, hitDiskIds[diskId], matID, myLocalData, pGlobalData_, rngState); } } else { // Triangle Geometry - single hit particle->surfaceCollision( - rayWeight, rayDir, geomNormal, rayHit.hit.primID, + rayWeight, rayDirection, geomNormal, rayHit.hit.primID, geometry_.getMaterialId(rayHit.hit.primID), myLocalData, pGlobalData_, rngState); } - // get sticking probability and reflection direction + // get sticking probability (1) and reflected direction (2) const auto stickingDirection = particle->surfaceReflection( - rayWeight, rayDir, geomNormal, rayHit.hit.primID, + rayWeight, rayDirection, geomNormal, rayHit.hit.primID, geometry_.getMaterialId(rayHit.hit.primID), pGlobalData_, rngState); @@ -331,7 +324,7 @@ template class TraceKernel { // Update ray direction and origin fillRayPosition(rayHit.ray, hitPoint); - fillRayDirection(rayHit.ray, stickingDirection.second); + fillRayDirection(rayHit.ray, stickingDirection.second); } while (reflect); totalBoundaryHits += boundaryHits; diff --git a/include/viennaray/rayUtil.hpp b/include/viennaray/rayUtil.hpp index 08f7f9f..1973fb9 100644 --- a/include/viennaray/rayUtil.hpp +++ b/include/viennaray/rayUtil.hpp @@ -201,61 +201,62 @@ getTraceSettings(TraceDirection sourceDir) { return set; } -template -void fillRayDirection(RTCRay &ray, const Vec3D &direction, - const float time = 0.0f) { +template +inline void fillRayDirection(RTCRay &ray, Vec3D direction, + float time = 0.0f) { + static_assert(std::is_same_v || std::is_same_v, + "fillRayDirection: T must be float or double"); + + if constexpr (D == 2) { + if (direction[2] != T(0)) { + direction[2] = T(0); + Normalize(direction); + } + } + #ifdef ARCH_X86 reinterpret_cast<__m128 &>(ray.dir_x) = _mm_set_ps( time, static_cast(direction[2]), static_cast(direction[1]), static_cast(direction[0])); #else - ray.dir_x = (float)direction[0]; - ray.dir_y = (float)direction[1]; - ray.dir_z = (float)direction[2]; + ray.dir_x = static_cast(direction[0]); + ray.dir_y = static_cast(direction[1]); + ray.dir_z = static_cast(direction[2]); ray.time = time; #endif } template -void fillRayPosition(RTCRay &ray, const Vec3D &origin, - const float tnear = 1e-4f) { +inline void fillRayPosition(RTCRay &ray, const Vec3D &origin, + float tnear = 1e-4f) { + static_assert(std::is_same_v || std::is_same_v, + "fillRayPosition: T must be float or double"); + #ifdef ARCH_X86 reinterpret_cast<__m128 &>(ray) = _mm_set_ps(tnear, static_cast(origin[2]), static_cast(origin[1]), static_cast(origin[0])); #else - ray.org_x = (float)origin[0]; - ray.org_y = (float)origin[1]; - ray.org_z = (float)origin[2]; + ray.org_x = static_cast(origin[0]); + ray.org_y = static_cast(origin[1]); + ray.org_z = static_cast(origin[2]); ray.tnear = tnear; #endif } -template <> -inline void fillRayDirection(RTCRay &ray, const Vec3D &direction, - const float time) { -#ifdef ARCH_X86 - reinterpret_cast<__m128 &>(ray.dir_x) = - _mm_set_ps(time, direction[2], direction[1], direction[0]); -#else - ray.dir_x = direction[0]; - ray.dir_y = direction[1]; - ray.dir_z = direction[2]; - ray.time = time; -#endif -} +template inline void sincos(double x, T &s, T &c) { + static_assert(std::is_same_v || std::is_same_v, + "sincos only supports float or double"); -template <> -inline void fillRayPosition(RTCRay &ray, const Vec3D &origin, - const float tnear) { -#ifdef ARCH_X86 - reinterpret_cast<__m128 &>(ray) = - _mm_set_ps(tnear, origin[2], origin[1], origin[0]); +#if defined(__GLIBC__) + if constexpr (std::is_same_v) { + sincosf(x, &s, &c); + } else { + ::sincos(x, &s, &c); + } #else - ray.org_x = (float)origin[0]; - ray.org_y = (float)origin[1]; - ray.org_z = (float)origin[2]; - ray.tnear = tnear; + s = std::sin(x); + c = std::cos(x); #endif } From 8ca1d68a6d0c08f90764c965ad4c50489cb358b3 Mon Sep 17 00:00:00 2001 From: reiter Date: Tue, 3 Feb 2026 16:15:03 +0100 Subject: [PATCH 02/10] Use 3D direction in 2D angle calculation (GPU) --- .gitignore | 3 +- CMakeLists.txt | 2 +- gpu/examples/trenchLines.cpp | 15 +++-- gpu/pipelines/GeneralPipelineDisk.cu | 46 +++++++++----- gpu/pipelines/GeneralPipelineLine.cu | 49 +++++++++----- gpu/pipelines/GeneralPipelineTriangle.cu | 39 ++++++++---- gpu/pipelines/Particle.cuh | 4 +- include/viennaray/gpu/raygBoundary.hpp | 16 ++--- include/viennaray/gpu/raygPerRayData.hpp | 6 +- include/viennaray/gpu/raygReflection.hpp | 23 +++---- include/viennaray/gpu/raygSource.hpp | 81 +++++++++++------------- 11 files changed, 162 insertions(+), 122 deletions(-) diff --git a/.gitignore b/.gitignore index 753d924..a3284d3 100644 --- a/.gitignore +++ b/.gitignore @@ -13,4 +13,5 @@ CMakeSettings.json docs/doxygen/*.txt !docs/doxygen/doxygen_template.txt *_old* -.cache* \ No newline at end of file +.cache* +__pycache__/ \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt index 441a642..edf3abb 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -99,7 +99,7 @@ include("cmake/cpm.cmake") CPMAddPackage( NAME ViennaCore - VERSION 1.9.2 + VERSION 1.9.5 GIT_REPOSITORY "https://github.com/ViennaTools/ViennaCore" OPTIONS "VIENNACORE_USE_GPU ${VIENNARAY_USE_GPU}") diff --git a/gpu/examples/trenchLines.cpp b/gpu/examples/trenchLines.cpp index c38adb5..5741945 100644 --- a/gpu/examples/trenchLines.cpp +++ b/gpu/examples/trenchLines.cpp @@ -31,7 +31,7 @@ int main() { gpu::Particle particle; particle.name = "Particle"; - particle.sticking = 1.f; + particle.sticking = 0.5f; particle.dataLabels = {"particleFlux"}; particle.materialSticking[7] = 1.f; particle.materialSticking[1] = .1f; @@ -40,15 +40,17 @@ int main() { std::vector cMap = { {0, gpu::CallableSlot::COLLISION, "__direct_callable__particleCollision"}, {0, gpu::CallableSlot::REFLECTION, - "__direct_callable__particleReflection"}}; + "__direct_callable__particleReflectionConstSticking"}, + }; gpu::TraceLine tracer(context); tracer.setGeometry(mesh); tracer.setMaterialIds(materialIds); tracer.setCallables("ViennaRayCallableWrapper", context->modulePath); tracer.setParticleCallableMap({pMap, cMap}); - tracer.setNumberOfRaysPerPoint(5000); tracer.insertNextParticle(particle); + tracer.setNumberOfRaysPerPoint(5000); + tracer.setMaxBoundaryHits(10); tracer.prepareParticlePrograms(); tracer.apply(); @@ -73,16 +75,17 @@ int main() { triangleTracer.setMaterialIds(materialIds); triangleTracer.setCallables("ViennaRayCallableWrapper", context->modulePath); triangleTracer.setParticleCallableMap({pMap, cMap}); - triangleTracer.setNumberOfRaysPerPoint(5000); triangleTracer.insertNextParticle(particle); + triangleTracer.setNumberOfRaysPerPoint(5000); + triangleTracer.setMaxBoundaryHits(10); triangleTracer.prepareParticlePrograms(); triangleTracer.apply(); triangleTracer.normalizeResults(); - flux = triangleTracer.getFlux(0, 0); + auto fluxTriangles = triangleTracer.getFlux(0, 0); rayInternal::writeVTP("trenchLines_triFlux.vtp", triMesh.nodes, - triMesh.triangles, flux); + triMesh.triangles, fluxTriangles); return 0; } diff --git a/gpu/pipelines/GeneralPipelineDisk.cu b/gpu/pipelines/GeneralPipelineDisk.cu index 776da5c..b485883 100644 --- a/gpu/pipelines/GeneralPipelineDisk.cu +++ b/gpu/pipelines/GeneralPipelineDisk.cu @@ -32,11 +32,11 @@ extern "C" __global__ void __intersection__() { const float radius = sbtData->radius; bool valid = true; - float prodOfDirections = DotProduct(normal, prd->dir); + float prodOfDirections = DotProduct(normal, prd->traceDir); // Backface hits have to be reported so CH can let the ray through or kill the // ray if needed - // valid &= DotProduct(prd->dir, normal) <= 0.0f; + // valid &= DotProduct(prd->traceDir, normal) <= 0.0f; // Check if ray is not parallel to the plane valid &= fabsf(prodOfDirections) >= 1e-6f; @@ -45,7 +45,7 @@ extern "C" __global__ void __intersection__() { float t = (ddneg - DotProduct(normal, prd->pos)) / prodOfDirections; // Avoid negative t or self intersections valid &= t > optixGetRayTmin(); - const Vec3Df intersection = prd->pos + prd->dir * t; + const Vec3Df intersection = prd->pos + prd->traceDir * t; // Check if within disk radius const Vec3Df diff = intersection - diskOrigin; @@ -77,14 +77,14 @@ extern "C" __global__ void __closesthit__() { const Vec3Df &normal = sbtData->base.normal[primID]; // If closest hit was on backside, let it through once - if (DotProduct(prd->dir, normal) > 0.0f) { + if (DotProduct(prd->traceDir, normal) > 0.0f) { // If back was hit a second time, kill the ray if (prd->hitFromBack) { prd->rayWeight = 0.f; return; } prd->hitFromBack = true; - prd->pos = prd->pos + prd->tMin * prd->dir; + prd->pos = prd->pos + prd->tMin * prd->traceDir; return; } @@ -139,6 +139,8 @@ extern "C" __global__ void __closesthit__() { callableIndex(launchParams.particleType, CallableSlot::REFLECTION); optixDirectCall(callIdx, sbtData, prd); + + prd->numReflections++; } } @@ -153,15 +155,22 @@ extern "C" __global__ void __raygen__() { // per-ray data PerRayData prd; // each ray has its own RNG state - initializeRNGState(&prd, linearLaunchIndex, launchParams.seed); + initializeRNGState(prd, linearLaunchIndex, launchParams.seed); // initialize ray position and direction - initializeRayPosition(&prd, launchParams.source, launchParams.D); + initializeRayPosition(prd, launchParams.source, launchParams.D); if (launchParams.source.customDirectionBasis) { - initializeRayDirection(&prd, launchParams.cosineExponent, - launchParams.source.directionBasis, launchParams.D); + initializeRayDirection(prd, launchParams.cosineExponent, + launchParams.source.directionBasis); + prd.traceDir = prd.dir; } else { - initializeRayDirection(&prd, launchParams.cosineExponent, launchParams.D); + initializeRayDirection(prd, launchParams.cosineExponent); + prd.traceDir = prd.dir; + if (launchParams.D == 2) { + // fold z into y for 2D + prd.dir[1] = prd.dir[2]; + prd.traceDir[1] = prd.traceDir[2]; + } } unsigned callIdx = @@ -175,12 +184,17 @@ extern "C" __global__ void __raygen__() { unsigned int hintBitLength = 2; while (continueRay(launchParams, prd)) { + if (launchParams.D == 2) { + prd.traceDir[2] = 0.f; + viennacore::Normalize(prd.traceDir); + } optixTraverse(launchParams.traversable, // traversable GAS make_float3(prd.pos[0], prd.pos[1], prd.pos[2]), // origin - make_float3(prd.dir[0], prd.dir[1], prd.dir[2]), // direction - 1e-4f, // tmin - 1e20f, // tmax - 0.0f, // rayTime + make_float3(prd.traceDir[0], prd.traceDir[1], + prd.traceDir[2]), // direction + 1e-4f, // tmin + 1e20f, // tmax + 0.0f, // rayTime OptixVisibilityMask(255), OPTIX_RAY_FLAG_DISABLE_ANYHIT, // OPTIX_RAY_FLAG_NONE, 0, // SBT offset @@ -198,7 +212,7 @@ extern "C" __global__ void __raygen__() { } optixReorder(hint, hintBitLength); optixInvoke(u0, u1); - prd.totalCount = 0; // Reset PerRayData - prd.numReflections++; + prd.totalCount = 0; // Reset PerRayData + prd.traceDir = prd.dir; // Update traceDir for the next iteration } } diff --git a/gpu/pipelines/GeneralPipelineLine.cu b/gpu/pipelines/GeneralPipelineLine.cu index 0e0d258..47a023b 100644 --- a/gpu/pipelines/GeneralPipelineLine.cu +++ b/gpu/pipelines/GeneralPipelineLine.cu @@ -16,7 +16,7 @@ using namespace viennaray::gpu; -extern "C" __constant__ viennaray::gpu::LaunchParams launchParams; +extern "C" __constant__ LaunchParams launchParams; extern "C" __global__ void __intersection__() { const HitSBTDataLine *sbtData = @@ -32,7 +32,8 @@ extern "C" __global__ void __intersection__() { const Vec3Df &p1 = sbtData->nodes[idx[1]]; Vec3Df lineDir = p1 - p0; - float d = 1.f / (prd->dir[0] * lineDir[1] - prd->dir[1] * lineDir[0]); + float d = + 1.f / (prd->traceDir[0] * lineDir[1] - prd->traceDir[1] * lineDir[0]); bool valid = true; @@ -41,9 +42,9 @@ extern "C" __global__ void __intersection__() { // float t = CrossProduct((p0 - prd->pos), lineDir)[2] * d; valid &= t > 1e-5f; - float s = - d * (p0ToRayOrigin[0] * prd->dir[1] - p0ToRayOrigin[1] * prd->dir[0]); - // float s = CrossProduct((p0 - prd->pos), prd->dir)[2] * d; + float s = d * (p0ToRayOrigin[0] * prd->traceDir[1] - + p0ToRayOrigin[1] * prd->traceDir[0]); + // float s = CrossProduct((p0 - prd->pos), prd->traceDir)[2] * d; valid &= s > 1e-5f && s < 1.0f - 1e-5f; if (valid) { @@ -81,6 +82,8 @@ extern "C" __global__ void __closesthit__() { callableIndex(launchParams.particleType, CallableSlot::REFLECTION); optixDirectCall(callIdx, sbtData, prd); + + prd->numReflections++; } } @@ -95,15 +98,20 @@ extern "C" __global__ void __raygen__() { // per-ray data PerRayData prd; // each ray has its own RNG state - initializeRNGState(&prd, linearLaunchIndex, launchParams.seed); + initializeRNGState(prd, linearLaunchIndex, launchParams.seed); // initialize ray position and direction - initializeRayPosition(&prd, launchParams.source, launchParams.D); + initializeRayPosition(prd, launchParams.source, launchParams.D); if (launchParams.source.customDirectionBasis) { - initializeRayDirection(&prd, launchParams.cosineExponent, - launchParams.source.directionBasis, launchParams.D); + initializeRayDirection(prd, launchParams.cosineExponent, + launchParams.source.directionBasis); } else { - initializeRayDirection(&prd, launchParams.cosineExponent, launchParams.D); + initializeRayDirection(prd, launchParams.cosineExponent); + if (launchParams.D == 2) { + // fold z into y for 2D + prd.dir[1] = prd.dir[2]; + prd.traceDir[1] = prd.traceDir[2]; + } } unsigned callIdx = @@ -117,12 +125,22 @@ extern "C" __global__ void __raygen__() { unsigned int hintBitLength = 2; while (continueRay(launchParams, prd)) { + if (launchParams.D == 2) { + prd.traceDir[2] = 0.f; + viennacore::Normalize(prd.traceDir); + } + // printf("Idx: %u, pos: (%f, %f, %f), dir: (%f, %f, %f), traceDir: " + // "(%f, %f, %f), r: %u\n", + // linearLaunchIndex, prd.pos[0], prd.pos[1], prd.pos[2], prd.dir[0], + // prd.dir[1], prd.dir[2], prd.traceDir[0], prd.traceDir[1], + // prd.traceDir[2], prd.numReflections); optixTraverse(launchParams.traversable, // traversable GAS make_float3(prd.pos[0], prd.pos[1], prd.pos[2]), // origin - make_float3(prd.dir[0], prd.dir[1], prd.dir[2]), // direction - 1e-4f, // tmin - 1e20f, // tmax - 0.0f, // rayTime + make_float3(prd.traceDir[0], prd.traceDir[1], + prd.traceDir[2]), // direction + 1e-4f, // tmin + 1e20f, // tmax + 0.0f, // rayTime OptixVisibilityMask(255), OPTIX_RAY_FLAG_DISABLE_ANYHIT, // OPTIX_RAY_FLAG_NONE, 0, // SBT offset @@ -140,7 +158,6 @@ extern "C" __global__ void __raygen__() { } optixReorder(hint, hintBitLength); optixInvoke(u0, u1); - prd.totalCount = 0; // Reset PerRayData - prd.numReflections++; + prd.traceDir = prd.dir; // Update traceDir for the next iteration } } diff --git a/gpu/pipelines/GeneralPipelineTriangle.cu b/gpu/pipelines/GeneralPipelineTriangle.cu index ebee29d..da926eb 100644 --- a/gpu/pipelines/GeneralPipelineTriangle.cu +++ b/gpu/pipelines/GeneralPipelineTriangle.cu @@ -18,7 +18,7 @@ using namespace viennaray::gpu; -extern "C" __constant__ viennaray::gpu::LaunchParams launchParams; +extern "C" __constant__ LaunchParams launchParams; extern "C" __global__ void __closesthit__() { const HitSBTDataTriangle *sbtData = @@ -51,6 +51,7 @@ extern "C" __global__ void __closesthit__() { callableIndex(launchParams.particleType, CallableSlot::REFLECTION); optixDirectCall( callIdx, sbtData, prd); + prd->numReflections++; } } @@ -65,15 +66,20 @@ extern "C" __global__ void __raygen__() { // per-ray data PerRayData prd; // each ray has its own RNG state - initializeRNGState(&prd, linearLaunchIndex, launchParams.seed); + initializeRNGState(prd, linearLaunchIndex, launchParams.seed); // initialize ray position and direction - initializeRayPosition(&prd, launchParams.source, launchParams.D); + initializeRayPosition(prd, launchParams.source, launchParams.D); if (launchParams.source.customDirectionBasis) { - initializeRayDirection(&prd, launchParams.cosineExponent, - launchParams.source.directionBasis, launchParams.D); + initializeRayDirection(prd, launchParams.cosineExponent, + launchParams.source.directionBasis); } else { - initializeRayDirection(&prd, launchParams.cosineExponent, launchParams.D); + initializeRayDirection(prd, launchParams.cosineExponent); + if (launchParams.D == 2) { + // fold z into y for 2D + prd.dir[1] = prd.dir[2]; + prd.traceDir[1] = prd.traceDir[2]; + } } unsigned callIdx = @@ -87,12 +93,22 @@ extern "C" __global__ void __raygen__() { unsigned int hintBitLength = 2; while (continueRay(launchParams, prd)) { + if (launchParams.D == 2) { + prd.traceDir[2] = 0.f; + viennacore::Normalize(prd.traceDir); + } + // printf("Idx: %u, pos: (%f, %f, %f), dir: (%f, %f, %f), traceDir: " + // "(%f, %f, %f), r: %u\n", + // linearLaunchIndex, prd.pos[0], prd.pos[1], prd.pos[2], prd.dir[0], + // prd.dir[1], prd.dir[2], prd.traceDir[0], prd.traceDir[1], + // prd.traceDir[2], prd.numReflections); optixTraverse(launchParams.traversable, // traversable GAS make_float3(prd.pos[0], prd.pos[1], prd.pos[2]), // origin - make_float3(prd.dir[0], prd.dir[1], prd.dir[2]), // direction - 1e-4f, // tmin - 1e20f, // tmax - 0.0f, // rayTime + make_float3(prd.traceDir[0], prd.traceDir[1], + prd.traceDir[2]), // direction + 1e-4f, // tmin + 1e20f, // tmax + 0.0f, // rayTime OptixVisibilityMask(255), OPTIX_RAY_FLAG_DISABLE_ANYHIT, // OPTIX_RAY_FLAG_NONE, 0, // SBT offset @@ -111,7 +127,8 @@ extern "C" __global__ void __raygen__() { } optixReorder(hint, hintBitLength); optixInvoke(u0, u1); - prd.numReflections++; + prd.traceDir = prd.dir; // Update traceDir for the next iteration + #ifdef COUNT_RAYS int *counter = reinterpret_cast(launchParams.customData); atomicAdd(counter, 1); diff --git a/gpu/pipelines/Particle.cuh b/gpu/pipelines/Particle.cuh index 7af9e64..344cdaf 100644 --- a/gpu/pipelines/Particle.cuh +++ b/gpu/pipelines/Particle.cuh @@ -27,7 +27,7 @@ particleReflection(const void *sbtData, viennaray::gpu::PerRayData *prd) { int materialId = launchParams.materialIds[prd->primID]; prd->rayWeight -= prd->rayWeight * launchParams.materialSticking[materialId]; auto geoNormal = viennaray::gpu::computeNormal(sbtData, prd->primID); - viennaray::gpu::diffuseReflection(prd, geoNormal, launchParams.D); + viennaray::gpu::diffuseReflection(prd, geoNormal); } __forceinline__ __device__ void @@ -35,7 +35,7 @@ particleReflectionConstSticking(const void *sbtData, viennaray::gpu::PerRayData *prd) { prd->rayWeight -= prd->rayWeight * launchParams.sticking; auto geoNormal = viennaray::gpu::computeNormal(sbtData, prd->primID); - viennaray::gpu::diffuseReflection(prd, geoNormal, launchParams.D); + viennaray::gpu::diffuseReflection(prd, geoNormal); } __forceinline__ __device__ void particleInit(viennaray::gpu::PerRayData *prd) { diff --git a/include/viennaray/gpu/raygBoundary.hpp b/include/viennaray/gpu/raygBoundary.hpp index 305e41b..9980b82 100644 --- a/include/viennaray/gpu/raygBoundary.hpp +++ b/include/viennaray/gpu/raygBoundary.hpp @@ -22,20 +22,20 @@ reflectFromBoundary(PerRayData *prd, const SBTData *hsd, const int D) { prd->numBoundaryHits++; if constexpr (std::is_same_v) { - prd->pos = - prd->pos + prd->dir * (optixGetRayTmax() - launchParams.tThreshold); + prd->pos = prd->pos + + prd->traceDir * (optixGetRayTmax() - launchParams.tThreshold); if (primID == 0 || primID == 1) { prd->dir[0] -= 2 * prd->dir[0]; // x boundary } else if ((primID == 2 || primID == 3) && D == 3) { prd->dir[1] -= 2 * prd->dir[1]; // y boundary } } else if constexpr (std::is_same_v) { - prd->pos = prd->pos + prd->dir * optixGetRayTmax(); + prd->pos = prd->pos + prd->traceDir * optixGetRayTmax(); unsigned dim = primID / 4; prd->dir[dim] -= 2 * prd->dir[dim]; prd->pos[dim] = hsd->vertex[hsd->index[primID][0]][dim]; } else if constexpr (std::is_same_v) { - prd->pos = prd->pos + prd->dir * optixGetRayTmax(); + prd->pos = prd->pos + prd->traceDir * optixGetRayTmax(); if (primID == 0 || primID == 1) // x boundary prd->dir[0] -= 2 * prd->dir[0]; } @@ -49,8 +49,8 @@ applyPeriodicBoundary(PerRayData *prd, const SBTData *hsd, const int D) { prd->numBoundaryHits++; if constexpr (std::is_same_v) { - prd->pos = - prd->pos + prd->dir * (optixGetRayTmax() - launchParams.tThreshold); + prd->pos = prd->pos + + prd->traceDir * (optixGetRayTmax() - launchParams.tThreshold); if (primID == 0) { // xmin prd->pos[0] = hsd->point[1][0]; } else if (primID == 1) { // xmax @@ -61,11 +61,11 @@ applyPeriodicBoundary(PerRayData *prd, const SBTData *hsd, const int D) { prd->pos[1] = hsd->point[2][1]; } } else if constexpr (std::is_same_v) { - prd->pos = prd->pos + prd->dir * optixGetRayTmax(); + prd->pos = prd->pos + prd->traceDir * optixGetRayTmax(); unsigned dim = primID / 4; prd->pos[dim] = hsd->vertex[hsd->index[primID ^ 2][0]][dim]; } else if constexpr (std::is_same_v) { - prd->pos = prd->pos + prd->dir * optixGetRayTmax(); + prd->pos = prd->pos + prd->traceDir * optixGetRayTmax(); if (primID == 0) { // xmin prd->pos[0] = hsd->nodes[1][0]; } else if (primID == 1) { // xmax diff --git a/include/viennaray/gpu/raygPerRayData.hpp b/include/viennaray/gpu/raygPerRayData.hpp index 45e73da..afdbe42 100644 --- a/include/viennaray/gpu/raygPerRayData.hpp +++ b/include/viennaray/gpu/raygPerRayData.hpp @@ -17,6 +17,8 @@ struct PerRayData { // Position and direction Vec3Df pos; Vec3Df dir; + Vec3Df traceDir; // direction for calculating the intersection, which can be + // different in 2D from the actual ray direction // Simulation specific data float rayWeight = 1.f; @@ -65,10 +67,10 @@ static __device__ __forceinline__ PerRayData *getPRD() { } static __device__ __forceinline__ void -initializeRNGState(PerRayData *prd, unsigned int linearLaunchIndex, +initializeRNGState(PerRayData &prd, unsigned int linearLaunchIndex, unsigned int seed) { auto rngSeed = tea<3>(linearLaunchIndex, seed); - curand_init(rngSeed, 0, 0, &prd->RNGstate); + curand_init(rngSeed, 0, 0, &prd.RNGstate); } #endif diff --git a/include/viennaray/gpu/raygReflection.hpp b/include/viennaray/gpu/raygReflection.hpp index 64fbcac..abb8a9f 100644 --- a/include/viennaray/gpu/raygReflection.hpp +++ b/include/viennaray/gpu/raygReflection.hpp @@ -60,7 +60,7 @@ __device__ __forceinline__ Vec3Df getNormal(const void *sbtData, static __device__ __forceinline__ void specularReflection(PerRayData *prd, const Vec3Df &geoNormal) { #ifndef VIENNARAY_TEST - prd->pos = prd->pos + prd->tMin * prd->dir; + prd->pos = prd->pos + prd->tMin * prd->traceDir; #endif prd->dir = prd->dir - (2 * DotProduct(prd->dir, geoNormal)) * geoNormal; } @@ -76,31 +76,27 @@ static __device__ Vec3Df PickRandomPointOnUnitSphere(CudaRNG *state) { return Vec3Df{r * c, r * s, z}; } -static __device__ void -diffuseReflection(PerRayData *prd, const Vec3Df &geoNormal, const uint8_t D) { +static __device__ void diffuseReflection(PerRayData *prd, + const Vec3Df &geoNormal) { #ifndef VIENNARAY_TEST - prd->pos = prd->pos + prd->tMin * prd->dir; + prd->pos = prd->pos + prd->tMin * prd->traceDir; #endif const Vec3Df randomDirection = PickRandomPointOnUnitSphere(&prd->RNGstate); prd->dir = geoNormal + randomDirection; - - if (D == 2) - prd->dir[2] = 0.f; - Normalize(prd->dir); } -static __device__ void diffuseReflection(PerRayData *prd, const uint8_t D) { +static __device__ void diffuseReflection(PerRayData *prd) { const HitSBTDataDisk *sbtData = (const HitSBTDataDisk *)optixGetSbtDataPointer(); const Vec3Df geoNormal = computeNormal(sbtData, optixGetPrimitiveIndex()); - diffuseReflection(prd, geoNormal, D); + diffuseReflection(prd, geoNormal); } static __device__ __forceinline__ void conedCosineReflection(PerRayData *prd, const Vec3Df &geomNormal, - const float maxConeAngle, const uint8_t D) { + const float maxConeAngle) { // Calculate specular direction specularReflection(prd, geomNormal); @@ -108,7 +104,7 @@ conedCosineReflection(PerRayData *prd, const Vec3Df &geomNormal, return; } if (maxConeAngle >= M_PI_2f) { - diffuseReflection(prd, geomNormal, D); + diffuseReflection(prd, geomNormal); return; } @@ -153,9 +149,6 @@ conedCosineReflection(PerRayData *prd, const Vec3Df &geomNormal, if (dp <= 0.f) dir = dir - 2.f * dp * geomNormal; - if (D == 2) - dir[2] = 0.f; - Normalize(dir); prd->dir = dir; } diff --git a/include/viennaray/gpu/raygSource.hpp b/include/viennaray/gpu/raygSource.hpp index deec5e0..8cc02a0 100644 --- a/include/viennaray/gpu/raygSource.hpp +++ b/include/viennaray/gpu/raygSource.hpp @@ -25,33 +25,27 @@ getOrthonormalBasis(const Vec3Df &n) { return {n, t, b2}; } -__device__ __forceinline__ void -initializeRayDirection(PerRayData *prd, const float power, const uint16_t D) { +__device__ __forceinline__ void initializeRayDirection(PerRayData &prd, + const float power) { // source direction - const float4 u = curand_uniform4(&prd->RNGstate); // (0,1] + const float4 u = curand_uniform4(&prd.RNGstate); // (0,1] const float cosTheta = powf(u.w, 1.f / (power + 1.f)); const float sinTheta = sqrtf(max(0.f, 1.f - cosTheta * cosTheta)); float sinPhi, cosPhi; __sincosf(2.f * M_PIf * u.x, &sinPhi, &cosPhi); - prd->dir[0] = cosPhi * sinTheta; - if (D == 2) { - prd->dir[1] = -cosTheta; - prd->dir[2] = 0.f; - Normalize(prd->dir); - } else { - prd->dir[1] = sinPhi * sinTheta; - prd->dir[2] = -cosTheta; - // already normalized - } + prd.dir[0] = cosPhi * sinTheta; + prd.dir[1] = sinPhi * sinTheta; + prd.dir[2] = -cosTheta; + prd.traceDir = prd.dir; } __device__ __forceinline__ void -initializeRayDirection(PerRayData *prd, const float power, - const std::array &basis, const uint16_t D) { +initializeRayDirection(PerRayData &prd, const float power, + const std::array &basis) { // source direction do { - const float4 u = curand_uniform4(&prd->RNGstate); // (0,1] + const float4 u = curand_uniform4(&prd.RNGstate); // (0,1] const float cosTheta = powf(u.w, 1.f / (power + 1.f)); const float sinTheta = sqrtf(max(0.f, 1.f - cosTheta * cosTheta)); float sinPhi, cosPhi; @@ -61,52 +55,51 @@ initializeRayDirection(PerRayData *prd, const float power, float ry = cosPhi * sinTheta; float rz = sinPhi * sinTheta; - prd->dir = basis[0] * rx + basis[1] * ry + basis[2] * rz; - } while (prd->dir[2] >= 0.f); - - if (D == 2) - prd->dir[2] = 0.f; + prd.dir = basis[0] * rx + basis[1] * ry + basis[2] * rz; + } while (prd.dir[2] >= 0.f); - Normalize(prd->dir); + viennacore::Normalize(prd.dir); + prd.traceDir = prd.dir; } __device__ __forceinline__ void -initializeRayPosition(PerRayData *prd, const LaunchParams::SourcePlane &source, +initializeRayPosition(PerRayData &prd, const LaunchParams::SourcePlane &source, const uint16_t D) { - const float4 u = curand_uniform4(&prd->RNGstate); // (0,1] - prd->pos[0] = + const float4 u = curand_uniform4(&prd.RNGstate); // (0,1] + prd.pos[0] = source.minPoint[0] + u.x * (source.maxPoint[0] - source.minPoint[0]); if (D == 2) { - prd->pos[1] = source.planeHeight; - prd->pos[2] = 0.f; + prd.pos[1] = source.planeHeight; + prd.pos[2] = 0.f; } else { - prd->pos[1] = + prd.pos[1] = source.minPoint[1] + u.y * (source.maxPoint[1] - source.minPoint[1]); - prd->pos[2] = source.planeHeight; + prd.pos[2] = source.planeHeight; } } // This is slightly faster because there is only one call to curand_uniform4 __device__ __forceinline__ void -initializeRayPositionAndDirection(PerRayData *prd, - const LaunchParams *launchParams) { - const float4 u = curand_uniform4(&prd->RNGstate); // (0,1] - prd->pos[0] = launchParams->source.minPoint[0] + - u.x * (launchParams->source.maxPoint[0] - - launchParams->source.minPoint[0]); - prd->pos[1] = launchParams->source.minPoint[1] + - u.y * (launchParams->source.maxPoint[1] - - launchParams->source.minPoint[1]); - prd->pos[2] = launchParams->source.planeHeight; - - const float tt = powf(u.w, 1.f / (launchParams->cosineExponent + 1.f)); +initializeRayPositionAndDirection(PerRayData &prd, + const LaunchParams &launchParams) { + const float4 u = curand_uniform4(&prd.RNGstate); // (0,1] + prd.pos[0] = + launchParams.source.minPoint[0] + + u.x * (launchParams.source.maxPoint[0] - launchParams.source.minPoint[0]); + prd.pos[1] = + launchParams.source.minPoint[1] + + u.y * (launchParams.source.maxPoint[1] - launchParams.source.minPoint[1]); + prd.pos[2] = launchParams.source.planeHeight; + + const float tt = powf(u.w, 1.f / (launchParams.cosineExponent + 1.f)); float s, c; __sincosf(2.f * M_PIf * u.z, &s, &c); const float sqrt1mtt = sqrtf(1.f - tt * tt); - prd->dir[0] = c * sqrt1mtt; - prd->dir[1] = s * sqrt1mtt; - prd->dir[2] = -tt; + prd.dir[0] = c * sqrt1mtt; + prd.dir[1] = s * sqrt1mtt; + prd.dir[2] = -tt; + prd.traceDir = prd.dir; } } // namespace viennaray::gpu #endif From 230eb73dcd14605312d8232562f8509a1b131f77 Mon Sep 17 00:00:00 2001 From: reiter Date: Wed, 4 Feb 2026 11:25:12 +0100 Subject: [PATCH 03/10] Use general ray init function --- gpu/examples/trenchTriangles.cpp | 7 ++++ gpu/pipelines/GeneralPipelineDisk.cu | 25 ++------------ gpu/pipelines/GeneralPipelineLine.cu | 28 ++------------- gpu/pipelines/GeneralPipelineTriangle.cu | 30 ++-------------- include/viennaray/gpu/raygSource.hpp | 44 +++++++++++++++--------- 5 files changed, 40 insertions(+), 94 deletions(-) diff --git a/gpu/examples/trenchTriangles.cpp b/gpu/examples/trenchTriangles.cpp index 740f150..43d75cd 100644 --- a/gpu/examples/trenchTriangles.cpp +++ b/gpu/examples/trenchTriangles.cpp @@ -57,10 +57,17 @@ int main(int argc, char **argv) { tracer.setParameters(rayCountBuffer.dPointer()); #endif + Timer timer; + timer.start(); tracer.apply(); tracer.normalizeResults(); + timer.finish(); + auto flux = tracer.getFlux(0, 0); + std::cout << "Tracing time: " << timer.currentDuration / 1e9 << " seconds." + << std::endl; + rayInternal::writeVTP( "trenchTriangles_triMesh.vtp", mesh.nodes, mesh.triangles, flux); diff --git a/gpu/pipelines/GeneralPipelineDisk.cu b/gpu/pipelines/GeneralPipelineDisk.cu index b485883..6349989 100644 --- a/gpu/pipelines/GeneralPipelineDisk.cu +++ b/gpu/pipelines/GeneralPipelineDisk.cu @@ -158,20 +158,7 @@ extern "C" __global__ void __raygen__() { initializeRNGState(prd, linearLaunchIndex, launchParams.seed); // initialize ray position and direction - initializeRayPosition(prd, launchParams.source, launchParams.D); - if (launchParams.source.customDirectionBasis) { - initializeRayDirection(prd, launchParams.cosineExponent, - launchParams.source.directionBasis); - prd.traceDir = prd.dir; - } else { - initializeRayDirection(prd, launchParams.cosineExponent); - prd.traceDir = prd.dir; - if (launchParams.D == 2) { - // fold z into y for 2D - prd.dir[1] = prd.dir[2]; - prd.traceDir[1] = prd.traceDir[2]; - } - } + initializeRayPositionAndDirection(prd, launchParams); unsigned callIdx = callableIndex(launchParams.particleType, CallableSlot::INIT); @@ -201,15 +188,7 @@ extern "C" __global__ void __raygen__() { 1, // SBT stride 0, // missSBTIndex u0, u1); // Payload - unsigned int hint = 0; - if (prd.rayWeight < launchParams.rayWeightThreshold || prd.energy < 0.f) { - hint |= (1 << 0); - } - if (optixHitObjectIsHit()) { - const HitSBTDataDisk *hitData = reinterpret_cast( - optixHitObjectGetSbtDataPointer()); - hint |= hitData->base.isBoundary << 1; - } + unsigned int hint = getCoherenceHint(prd, launchParams); optixReorder(hint, hintBitLength); optixInvoke(u0, u1); prd.totalCount = 0; // Reset PerRayData diff --git a/gpu/pipelines/GeneralPipelineLine.cu b/gpu/pipelines/GeneralPipelineLine.cu index 47a023b..1f929d1 100644 --- a/gpu/pipelines/GeneralPipelineLine.cu +++ b/gpu/pipelines/GeneralPipelineLine.cu @@ -101,18 +101,7 @@ extern "C" __global__ void __raygen__() { initializeRNGState(prd, linearLaunchIndex, launchParams.seed); // initialize ray position and direction - initializeRayPosition(prd, launchParams.source, launchParams.D); - if (launchParams.source.customDirectionBasis) { - initializeRayDirection(prd, launchParams.cosineExponent, - launchParams.source.directionBasis); - } else { - initializeRayDirection(prd, launchParams.cosineExponent); - if (launchParams.D == 2) { - // fold z into y for 2D - prd.dir[1] = prd.dir[2]; - prd.traceDir[1] = prd.traceDir[2]; - } - } + initializeRayPositionAndDirection(prd, launchParams); unsigned callIdx = callableIndex(launchParams.particleType, CallableSlot::INIT); @@ -129,11 +118,6 @@ extern "C" __global__ void __raygen__() { prd.traceDir[2] = 0.f; viennacore::Normalize(prd.traceDir); } - // printf("Idx: %u, pos: (%f, %f, %f), dir: (%f, %f, %f), traceDir: " - // "(%f, %f, %f), r: %u\n", - // linearLaunchIndex, prd.pos[0], prd.pos[1], prd.pos[2], prd.dir[0], - // prd.dir[1], prd.dir[2], prd.traceDir[0], prd.traceDir[1], - // prd.traceDir[2], prd.numReflections); optixTraverse(launchParams.traversable, // traversable GAS make_float3(prd.pos[0], prd.pos[1], prd.pos[2]), // origin make_float3(prd.traceDir[0], prd.traceDir[1], @@ -147,15 +131,7 @@ extern "C" __global__ void __raygen__() { 1, // SBT stride 0, // missSBTIndex u0, u1); // Payload - unsigned int hint = 0; - if (prd.rayWeight < launchParams.rayWeightThreshold || prd.energy < 0.f) { - hint |= (1 << 0); - } - if (optixHitObjectIsHit()) { - const HitSBTDataLine *hitData = reinterpret_cast( - optixHitObjectGetSbtDataPointer()); - hint |= hitData->base.isBoundary << 1; - } + unsigned int hint = getCoherenceHint(prd, launchParams); optixReorder(hint, hintBitLength); optixInvoke(u0, u1); prd.traceDir = prd.dir; // Update traceDir for the next iteration diff --git a/gpu/pipelines/GeneralPipelineTriangle.cu b/gpu/pipelines/GeneralPipelineTriangle.cu index da926eb..7fdf6ad 100644 --- a/gpu/pipelines/GeneralPipelineTriangle.cu +++ b/gpu/pipelines/GeneralPipelineTriangle.cu @@ -69,18 +69,7 @@ extern "C" __global__ void __raygen__() { initializeRNGState(prd, linearLaunchIndex, launchParams.seed); // initialize ray position and direction - initializeRayPosition(prd, launchParams.source, launchParams.D); - if (launchParams.source.customDirectionBasis) { - initializeRayDirection(prd, launchParams.cosineExponent, - launchParams.source.directionBasis); - } else { - initializeRayDirection(prd, launchParams.cosineExponent); - if (launchParams.D == 2) { - // fold z into y for 2D - prd.dir[1] = prd.dir[2]; - prd.traceDir[1] = prd.traceDir[2]; - } - } + initializeRayPositionAndDirection(prd, launchParams); unsigned callIdx = callableIndex(launchParams.particleType, CallableSlot::INIT); @@ -97,11 +86,6 @@ extern "C" __global__ void __raygen__() { prd.traceDir[2] = 0.f; viennacore::Normalize(prd.traceDir); } - // printf("Idx: %u, pos: (%f, %f, %f), dir: (%f, %f, %f), traceDir: " - // "(%f, %f, %f), r: %u\n", - // linearLaunchIndex, prd.pos[0], prd.pos[1], prd.pos[2], prd.dir[0], - // prd.dir[1], prd.dir[2], prd.traceDir[0], prd.traceDir[1], - // prd.traceDir[2], prd.numReflections); optixTraverse(launchParams.traversable, // traversable GAS make_float3(prd.pos[0], prd.pos[1], prd.pos[2]), // origin make_float3(prd.traceDir[0], prd.traceDir[1], @@ -115,20 +99,10 @@ extern "C" __global__ void __raygen__() { 1, // SBT stride 0, // missSBTIndex u0, u1); // Payload - unsigned int hint = 0; - if (prd.rayWeight < launchParams.rayWeightThreshold || prd.energy < 0.f) { - hint |= (1 << 0); - } - if (optixHitObjectIsHit()) { - const HitSBTDataTriangle *hitData = - reinterpret_cast( - optixHitObjectGetSbtDataPointer()); - hint |= hitData->base.isBoundary << 1; - } + unsigned int hint = getCoherenceHint(prd, launchParams); optixReorder(hint, hintBitLength); optixInvoke(u0, u1); prd.traceDir = prd.dir; // Update traceDir for the next iteration - #ifdef COUNT_RAYS int *counter = reinterpret_cast(launchParams.customData); atomicAdd(counter, 1); diff --git a/include/viennaray/gpu/raygSource.hpp b/include/viennaray/gpu/raygSource.hpp index 8cc02a0..dd3aede 100644 --- a/include/viennaray/gpu/raygSource.hpp +++ b/include/viennaray/gpu/raygSource.hpp @@ -83,23 +83,33 @@ initializeRayPosition(PerRayData &prd, const LaunchParams::SourcePlane &source, __device__ __forceinline__ void initializeRayPositionAndDirection(PerRayData &prd, const LaunchParams &launchParams) { - const float4 u = curand_uniform4(&prd.RNGstate); // (0,1] - prd.pos[0] = - launchParams.source.minPoint[0] + - u.x * (launchParams.source.maxPoint[0] - launchParams.source.minPoint[0]); - prd.pos[1] = - launchParams.source.minPoint[1] + - u.y * (launchParams.source.maxPoint[1] - launchParams.source.minPoint[1]); - prd.pos[2] = launchParams.source.planeHeight; - - const float tt = powf(u.w, 1.f / (launchParams.cosineExponent + 1.f)); - float s, c; - __sincosf(2.f * M_PIf * u.z, &s, &c); - const float sqrt1mtt = sqrtf(1.f - tt * tt); - prd.dir[0] = c * sqrt1mtt; - prd.dir[1] = s * sqrt1mtt; - prd.dir[2] = -tt; - prd.traceDir = prd.dir; + initializeRayPosition(prd, launchParams.source, launchParams.D); + if (launchParams.source.customDirectionBasis) { + initializeRayDirection(prd, launchParams.cosineExponent, + launchParams.source.directionBasis); + } else { + initializeRayDirection(prd, launchParams.cosineExponent); + if (launchParams.D == 2) { + // fold z into y for 2D + prd.dir[1] = prd.dir[2]; + prd.traceDir[1] = prd.traceDir[2]; + } + } } + +__device__ __forceinline__ unsigned +getCoherenceHint(PerRayData &prd, const LaunchParams &launchParams) { + unsigned int hint = 0; + if (prd.rayWeight < launchParams.rayWeightThreshold || prd.energy < 0.f) { + hint |= (1 << 0); + } + if (optixHitObjectIsHit()) { + const HitSBTDataDisk *hitData = reinterpret_cast( + optixHitObjectGetSbtDataPointer()); + hint |= hitData->base.isBoundary << 1; + } + return hint; +} + } // namespace viennaray::gpu #endif From d369243332c6283d3c5e30b7f663fb0c9627675f Mon Sep 17 00:00:00 2001 From: reiter Date: Wed, 4 Feb 2026 11:56:29 +0100 Subject: [PATCH 04/10] Seperate boundary hit shader --- gpu/pipelines/GeneralPipelineDisk.cu | 104 +++++++++++--------- gpu/pipelines/GeneralPipelineLine.cu | 51 +++++----- gpu/pipelines/GeneralPipelineTriangle.cu | 51 +++++----- include/viennaray/gpu/raygTrace.hpp | 16 ++- include/viennaray/gpu/raygTraceDisk.hpp | 3 +- include/viennaray/gpu/raygTraceLine.hpp | 3 +- include/viennaray/gpu/raygTraceTriangle.hpp | 3 +- 7 files changed, 136 insertions(+), 95 deletions(-) diff --git a/gpu/pipelines/GeneralPipelineDisk.cu b/gpu/pipelines/GeneralPipelineDisk.cu index 6349989..307e93c 100644 --- a/gpu/pipelines/GeneralPipelineDisk.cu +++ b/gpu/pipelines/GeneralPipelineDisk.cu @@ -88,59 +88,69 @@ extern "C" __global__ void __closesthit__() { return; } - if (sbtData->base.isBoundary) { - // This is effectively the miss shader - if (launchParams.D == 2 && - (primID == 2 || primID == 3)) { // bottom or top - ymin or ymax - prd->rayWeight = 0.0f; - return; + // ------------- NEIGHBOR FILTERING --------------- // + // Keep only hits close to tMin + prd->ISCount = 0; + for (int i = 0; i < prd->totalCount; ++i) { + if (fabsf(prd->tValues[i] - prd->tMin) < launchParams.tThreshold && + prd->ISCount < MAX_NEIGHBORS) { + prd->primIDs[prd->ISCount++] = prd->primIDs[i]; } - if (launchParams.D == 3 && - (primID == 4 || primID == 5)) { // bottom or top - zmin or zmax - prd->rayWeight = 0.0f; + } + + // ------------- SURFACE COLLISION --------------- // + unsigned callIdx = + callableIndex(launchParams.particleType, CallableSlot::COLLISION); + optixDirectCall(callIdx, sbtData, + prd); + + // ------------- REFLECTION --------------- // + callIdx = callableIndex(launchParams.particleType, CallableSlot::REFLECTION); + optixDirectCall(callIdx, sbtData, + prd); + + prd->numReflections++; +} + +extern "C" __global__ void __closesthit__boundary__() { + const HitSBTDataDisk *sbtData = + (const HitSBTDataDisk *)optixGetSbtDataPointer(); + PerRayData *prd = getPRD(); + + const unsigned int primID = optixGetPrimitiveIndex(); + prd->tMin = optixGetRayTmax() - launchParams.tThreshold; + prd->primID = primID; + + const Vec3Df &normal = sbtData->base.normal[primID]; + + // If closest hit was on backside, let it through once + if (DotProduct(prd->traceDir, normal) > 0.0f) { + // If back was hit a second time, kill the ray + if (prd->hitFromBack) { + prd->rayWeight = 0.f; return; } + prd->hitFromBack = true; + prd->pos = prd->pos + prd->tMin * prd->traceDir; + return; + } - if (launchParams.periodicBoundary) { - applyPeriodicBoundary(prd, sbtData, launchParams.D); - } else { - reflectFromBoundary(prd, sbtData, launchParams.D); - } + // This is effectively the miss shader + if (launchParams.D == 2 && + (primID == 2 || primID == 3)) { // bottom or top - ymin or ymax + prd->rayWeight = 0.0f; + return; + } + if (launchParams.D == 3 && + (primID == 4 || primID == 5)) { // bottom or top - zmin or zmax + prd->rayWeight = 0.0f; + return; + } + if (launchParams.periodicBoundary) { + applyPeriodicBoundary(prd, sbtData, launchParams.D); } else { - // ------------- NEIGHBOR FILTERING --------------- // - // Keep only hits close to tMin - prd->ISCount = 0; - for (int i = 0; i < prd->totalCount; ++i) { - if (fabsf(prd->tValues[i] - prd->tMin) < launchParams.tThreshold && - prd->ISCount < MAX_NEIGHBORS) { - prd->primIDs[prd->ISCount++] = prd->primIDs[i]; - } - } - - // // CPU like neighbor detection - // prd->ISCount = 0; - // for (int i = 0; i < prd->totalCount; ++i) { - // float distance = viennacore::Distance(sbtData->point[primID], - // sbtData->point[prd->primIDs[i]]); - // if (distance < 2 * sbtData->radius && prd->ISCount < MAX_NEIGHBORS) { - // prd->TIndex[prd->ISCount++] = prd->primIDs[i]; - // } - // } - - // ------------- SURFACE COLLISION --------------- // - unsigned callIdx = - callableIndex(launchParams.particleType, CallableSlot::COLLISION); - optixDirectCall(callIdx, - sbtData, prd); - - // ------------- REFLECTION --------------- // - callIdx = - callableIndex(launchParams.particleType, CallableSlot::REFLECTION); - optixDirectCall(callIdx, - sbtData, prd); - - prd->numReflections++; + reflectFromBoundary(prd, sbtData, launchParams.D); } } diff --git a/gpu/pipelines/GeneralPipelineLine.cu b/gpu/pipelines/GeneralPipelineLine.cu index 1f929d1..2914a8b 100644 --- a/gpu/pipelines/GeneralPipelineLine.cu +++ b/gpu/pipelines/GeneralPipelineLine.cu @@ -61,29 +61,36 @@ extern "C" __global__ void __closesthit__() { prd->tMin = optixGetRayTmax(); prd->primID = primID; - if (sbtData->base.isBoundary) { - if (launchParams.periodicBoundary) { - applyPeriodicBoundary(prd, sbtData, launchParams.D); - } else { - reflectFromBoundary(prd, sbtData, launchParams.D); - } + prd->ISCount = 1; + prd->primIDs[0] = primID; + + // ------------- SURFACE COLLISION --------------- // + unsigned callIdx = + callableIndex(launchParams.particleType, CallableSlot::COLLISION); + optixDirectCall(callIdx, sbtData, + prd); + + // ------------- REFLECTION --------------- // + callIdx = callableIndex(launchParams.particleType, CallableSlot::REFLECTION); + optixDirectCall(callIdx, sbtData, + prd); + + prd->numReflections++; +} + +extern "C" __global__ void __closesthit__boundary__() { + const HitSBTDataLine *sbtData = + (const HitSBTDataLine *)optixGetSbtDataPointer(); + PerRayData *prd = getPRD(); + + const unsigned int primID = optixGetPrimitiveIndex(); + prd->tMin = optixGetRayTmax(); + prd->primID = primID; + + if (launchParams.periodicBoundary) { + applyPeriodicBoundary(prd, sbtData, launchParams.D); } else { - prd->ISCount = 1; - prd->primIDs[0] = primID; - - // ------------- SURFACE COLLISION --------------- // - unsigned callIdx = - callableIndex(launchParams.particleType, CallableSlot::COLLISION); - optixDirectCall(callIdx, - sbtData, prd); - - // ------------- REFLECTION --------------- // - callIdx = - callableIndex(launchParams.particleType, CallableSlot::REFLECTION); - optixDirectCall(callIdx, - sbtData, prd); - - prd->numReflections++; + reflectFromBoundary(prd, sbtData, launchParams.D); } } diff --git a/gpu/pipelines/GeneralPipelineTriangle.cu b/gpu/pipelines/GeneralPipelineTriangle.cu index 7fdf6ad..8d3bb5c 100644 --- a/gpu/pipelines/GeneralPipelineTriangle.cu +++ b/gpu/pipelines/GeneralPipelineTriangle.cu @@ -29,29 +29,36 @@ extern "C" __global__ void __closesthit__() { prd->tMin = optixGetRayTmax(); prd->primID = primID; - if (sbtData->base.isBoundary) { - if (launchParams.periodicBoundary) { - applyPeriodicBoundary(prd, sbtData, launchParams.D); - } else { - reflectFromBoundary(prd, sbtData, launchParams.D); - } + prd->ISCount = 1; + prd->primIDs[0] = primID; + + // ------------- SURFACE COLLISION --------------- // + unsigned callIdx; + + callIdx = callableIndex(launchParams.particleType, CallableSlot::COLLISION); + optixDirectCall(callIdx, sbtData, prd); + + // ------------- REFLECTION --------------- // + callIdx = callableIndex(launchParams.particleType, CallableSlot::REFLECTION); + optixDirectCall(callIdx, + sbtData, prd); + prd->numReflections++; +} + +extern "C" __global__ void __closesthit__boundary__() { + const HitSBTDataTriangle *sbtData = + (const HitSBTDataTriangle *)optixGetSbtDataPointer(); + PerRayData *prd = getPRD(); + + const unsigned int primID = optixGetPrimitiveIndex(); + prd->tMin = optixGetRayTmax(); + prd->primID = primID; + + if (launchParams.periodicBoundary) { + applyPeriodicBoundary(prd, sbtData, launchParams.D); } else { - prd->ISCount = 1; - prd->primIDs[0] = primID; - - // ------------- SURFACE COLLISION --------------- // - unsigned callIdx; - - callIdx = callableIndex(launchParams.particleType, CallableSlot::COLLISION); - optixDirectCall(callIdx, sbtData, prd); - - // ------------- REFLECTION --------------- // - callIdx = - callableIndex(launchParams.particleType, CallableSlot::REFLECTION); - optixDirectCall( - callIdx, sbtData, prd); - prd->numReflections++; + reflectFromBoundary(prd, sbtData, launchParams.D); } } diff --git a/include/viennaray/gpu/raygTrace.hpp b/include/viennaray/gpu/raygTrace.hpp index c75b49a..240ea1b 100644 --- a/include/viennaray/gpu/raygTrace.hpp +++ b/include/viennaray/gpu/raygTrace.hpp @@ -631,8 +631,19 @@ template class Trace { pgDesc.hitgroup.moduleIS = module_; pgDesc.hitgroup.entryFunctionNameIS = entryFunctionNameIS.c_str(); } - createProgramGroup(&pgDesc, &pgOptions, &hitgroupPG_); + + std::string entryFunctionNameCHBound = "__closesthit__boundary__"; + OptixProgramGroupDesc pgDescBound = {}; + pgDescBound.kind = OPTIX_PROGRAM_GROUP_KIND_HITGROUP; + pgDescBound.hitgroup.moduleCH = module_; + pgDescBound.hitgroup.entryFunctionNameCH = entryFunctionNameCHBound.c_str(); + + if (geometryType_ != "Triangle") { + pgDescBound.hitgroup.moduleIS = module_; + pgDescBound.hitgroup.entryFunctionNameIS = entryFunctionNameIS.c_str(); + } + createProgramGroup(&pgDescBound, &pgOptions, &boundaryHitgroupPG_); } /// does all setup for the direct callables @@ -674,6 +685,8 @@ template class Trace { programGroups.push_back(raygenPG_); programGroups.push_back(missPG_); programGroups.push_back(hitgroupPG_); + if (boundaryHitgroupPG_) + programGroups.push_back(boundaryHitgroupPG_); for (auto const &directCallablePG : directCallablePGs_) { programGroups.push_back(directCallablePG); @@ -801,6 +814,7 @@ template class Trace { CudaBuffer missRecordBuffer_; OptixProgramGroup hitgroupPG_{}; CudaBuffer hitgroupRecordBuffer_; + OptixProgramGroup boundaryHitgroupPG_{}; std::vector directCallablePGs_; CudaBuffer directCallableRecordBuffer_; OptixShaderBindingTable shaderBindingTable_{}; diff --git a/include/viennaray/gpu/raygTraceDisk.hpp b/include/viennaray/gpu/raygTraceDisk.hpp index db442fb..1e9655f 100644 --- a/include/viennaray/gpu/raygTraceDisk.hpp +++ b/include/viennaray/gpu/raygTraceDisk.hpp @@ -199,7 +199,8 @@ template class TraceDisk final : public Trace { // boundary hitgroup if (!ignoreBoundary_) { HitgroupRecordDisk boundaryHitgroupRecord = {}; - optixSbtRecordPackHeader(this->hitgroupPG_, &boundaryHitgroupRecord); + optixSbtRecordPackHeader(this->boundaryHitgroupPG_, + &boundaryHitgroupRecord); boundaryHitgroupRecord.data.point = (Vec3Df *)diskGeometry_.boundaryPointBuffer.dPointer(); boundaryHitgroupRecord.data.base.normal = diff --git a/include/viennaray/gpu/raygTraceLine.hpp b/include/viennaray/gpu/raygTraceLine.hpp index 89f0444..b5d8b96 100644 --- a/include/viennaray/gpu/raygTraceLine.hpp +++ b/include/viennaray/gpu/raygTraceLine.hpp @@ -79,7 +79,8 @@ template class TraceLine final : public Trace { // boundary hitgroup if (!ignoreBoundary_) { HitgroupRecordLine boundaryHitgroupRecord = {}; - optixSbtRecordPackHeader(this->hitgroupPG_, &boundaryHitgroupRecord); + optixSbtRecordPackHeader(this->boundaryHitgroupPG_, + &boundaryHitgroupRecord); boundaryHitgroupRecord.data.nodes = (Vec3Df *)lineGeometry_.boundaryNodesBuffer.dPointer(); boundaryHitgroupRecord.data.lines = diff --git a/include/viennaray/gpu/raygTraceTriangle.hpp b/include/viennaray/gpu/raygTraceTriangle.hpp index 43219f2..2649a90 100644 --- a/include/viennaray/gpu/raygTraceTriangle.hpp +++ b/include/viennaray/gpu/raygTraceTriangle.hpp @@ -89,7 +89,8 @@ template class TraceTriangle final : public Trace { // boundary hitgroup if (!ignoreBoundary_) { HitgroupRecordTriangle boundaryHitgroupRecord = {}; - optixSbtRecordPackHeader(this->hitgroupPG_, &boundaryHitgroupRecord); + optixSbtRecordPackHeader(this->boundaryHitgroupPG_, + &boundaryHitgroupRecord); boundaryHitgroupRecord.data.vertex = (Vec3Df *)triangleGeometry_.boundaryVertexBuffer.dPointer(); boundaryHitgroupRecord.data.index = From 5124d47a6e4fc13ac7862abaa5cfc1fa923d6cc4 Mon Sep 17 00:00:00 2001 From: reiter Date: Wed, 4 Feb 2026 15:38:53 +0100 Subject: [PATCH 05/10] refactor: inline boundary functions --- gpu/pipelines/GeneralPipelineDisk.cu | 19 ++++++++----------- gpu/pipelines/GeneralPipelineLine.cu | 12 +++++------- gpu/pipelines/GeneralPipelineTriangle.cu | 16 +++++++--------- 3 files changed, 20 insertions(+), 27 deletions(-) diff --git a/gpu/pipelines/GeneralPipelineDisk.cu b/gpu/pipelines/GeneralPipelineDisk.cu index 307e93c..dfe106e 100644 --- a/gpu/pipelines/GeneralPipelineDisk.cu +++ b/gpu/pipelines/GeneralPipelineDisk.cu @@ -4,11 +4,9 @@ #define __CUDACC__ #endif -#include "raygBoundary.hpp" #include "raygCallableConfig.hpp" #include "raygLaunchParams.hpp" #include "raygPerRayData.hpp" -#include "raygReflection.hpp" #include "raygSBTRecords.hpp" #include "raygSource.hpp" @@ -123,14 +121,8 @@ extern "C" __global__ void __closesthit__boundary__() { const Vec3Df &normal = sbtData->base.normal[primID]; - // If closest hit was on backside, let it through once + // If closest hit was on backside of boundary, let it through if (DotProduct(prd->traceDir, normal) > 0.0f) { - // If back was hit a second time, kill the ray - if (prd->hitFromBack) { - prd->rayWeight = 0.f; - return; - } - prd->hitFromBack = true; prd->pos = prd->pos + prd->tMin * prd->traceDir; return; } @@ -147,10 +139,15 @@ extern "C" __global__ void __closesthit__boundary__() { return; } + // update ray position to hit point + prd->pos = prd->pos + prd->traceDir * prd->tMin; + + unsigned axis = primID / 2; if (launchParams.periodicBoundary) { - applyPeriodicBoundary(prd, sbtData, launchParams.D); + prd->pos[axis] = sbtData->point[primID ^ 1][axis]; // wrap to opposite side } else { - reflectFromBoundary(prd, sbtData, launchParams.D); + prd->dir[axis] -= 2 * prd->dir[axis]; // reflect + prd->pos[axis] = sbtData->point[primID][axis]; } } diff --git a/gpu/pipelines/GeneralPipelineLine.cu b/gpu/pipelines/GeneralPipelineLine.cu index 2914a8b..5cd60d0 100644 --- a/gpu/pipelines/GeneralPipelineLine.cu +++ b/gpu/pipelines/GeneralPipelineLine.cu @@ -4,11 +4,9 @@ #define __CUDACC__ #endif -#include "raygBoundary.hpp" #include "raygCallableConfig.hpp" #include "raygLaunchParams.hpp" #include "raygPerRayData.hpp" -#include "raygReflection.hpp" #include "raygSBTRecords.hpp" #include "raygSource.hpp" @@ -83,14 +81,14 @@ extern "C" __global__ void __closesthit__boundary__() { (const HitSBTDataLine *)optixGetSbtDataPointer(); PerRayData *prd = getPRD(); - const unsigned int primID = optixGetPrimitiveIndex(); - prd->tMin = optixGetRayTmax(); - prd->primID = primID; + // update ray position to hit point + prd->pos = prd->pos + prd->traceDir * optixGetRayTmax(); + const unsigned int primID = optixGetPrimitiveIndex(); if (launchParams.periodicBoundary) { - applyPeriodicBoundary(prd, sbtData, launchParams.D); + prd->pos[0] = sbtData->nodes[primID ^ 1][0]; // wrap around x-coordinate } else { - reflectFromBoundary(prd, sbtData, launchParams.D); + prd->dir[0] -= 2 * prd->dir[0]; // reflect } } diff --git a/gpu/pipelines/GeneralPipelineTriangle.cu b/gpu/pipelines/GeneralPipelineTriangle.cu index 8d3bb5c..bf0a1ef 100644 --- a/gpu/pipelines/GeneralPipelineTriangle.cu +++ b/gpu/pipelines/GeneralPipelineTriangle.cu @@ -4,11 +4,9 @@ #define __CUDACC__ #endif -#include "raygBoundary.hpp" #include "raygCallableConfig.hpp" #include "raygLaunchParams.hpp" #include "raygPerRayData.hpp" -#include "raygReflection.hpp" #include "raygSBTRecords.hpp" #include "raygSource.hpp" @@ -28,13 +26,11 @@ extern "C" __global__ void __closesthit__() { const unsigned int primID = optixGetPrimitiveIndex(); prd->tMin = optixGetRayTmax(); prd->primID = primID; - prd->ISCount = 1; prd->primIDs[0] = primID; // ------------- SURFACE COLLISION --------------- // unsigned callIdx; - callIdx = callableIndex(launchParams.particleType, CallableSlot::COLLISION); optixDirectCall(callIdx, sbtData, prd); @@ -51,14 +47,16 @@ extern "C" __global__ void __closesthit__boundary__() { (const HitSBTDataTriangle *)optixGetSbtDataPointer(); PerRayData *prd = getPRD(); - const unsigned int primID = optixGetPrimitiveIndex(); - prd->tMin = optixGetRayTmax(); - prd->primID = primID; + // update ray position to hit point + prd->pos = prd->pos + prd->traceDir * optixGetRayTmax(); + const unsigned int primID = optixGetPrimitiveIndex(); + unsigned dim = primID / 4; if (launchParams.periodicBoundary) { - applyPeriodicBoundary(prd, sbtData, launchParams.D); + prd->pos[dim] = sbtData->vertex[sbtData->index[primID ^ 2][0]][dim]; } else { - reflectFromBoundary(prd, sbtData, launchParams.D); + prd->dir[dim] -= 2 * prd->dir[dim]; + prd->pos[dim] = sbtData->vertex[sbtData->index[primID][0]][dim]; } } From b3e0ecf480e2f9afe8b3e048623d3ad8b8732a1d Mon Sep 17 00:00:00 2001 From: reiter Date: Wed, 4 Feb 2026 15:57:23 +0100 Subject: [PATCH 06/10] Small fixes --- gpu/pipelines/GeneralPipelineDisk.cu | 6 +----- gpu/pipelines/GeneralPipelineLine.cu | 6 ++---- gpu/pipelines/GeneralPipelineTriangle.cu | 10 ++-------- include/viennaray/gpu/raygSource.hpp | 7 ++++--- 4 files changed, 9 insertions(+), 20 deletions(-) diff --git a/gpu/pipelines/GeneralPipelineDisk.cu b/gpu/pipelines/GeneralPipelineDisk.cu index dfe106e..43b8d37 100644 --- a/gpu/pipelines/GeneralPipelineDisk.cu +++ b/gpu/pipelines/GeneralPipelineDisk.cu @@ -14,7 +14,7 @@ using namespace viennaray::gpu; -extern "C" __constant__ viennaray::gpu::LaunchParams launchParams; +extern "C" __constant__ LaunchParams launchParams; extern "C" __global__ void __intersection__() { const HitSBTDataDisk *sbtData = @@ -32,10 +32,6 @@ extern "C" __global__ void __intersection__() { bool valid = true; float prodOfDirections = DotProduct(normal, prd->traceDir); - // Backface hits have to be reported so CH can let the ray through or kill the - // ray if needed - // valid &= DotProduct(prd->traceDir, normal) <= 0.0f; - // Check if ray is not parallel to the plane valid &= fabsf(prodOfDirections) >= 1e-6f; diff --git a/gpu/pipelines/GeneralPipelineLine.cu b/gpu/pipelines/GeneralPipelineLine.cu index 5cd60d0..f21a144 100644 --- a/gpu/pipelines/GeneralPipelineLine.cu +++ b/gpu/pipelines/GeneralPipelineLine.cu @@ -22,7 +22,7 @@ extern "C" __global__ void __intersection__() { PerRayData *prd = getPRD(); // Get the index of the AABB box that was hit - const int primID = optixGetPrimitiveIndex(); + const unsigned int primID = optixGetPrimitiveIndex(); // Read geometric data from the primitive that is inside that AABB box const Vec2D &idx = (sbtData->lines)[primID]; @@ -37,12 +37,10 @@ extern "C" __global__ void __intersection__() { const Vec3Df p0ToRayOrigin = p0 - prd->pos; float t = d * (p0ToRayOrigin[0] * lineDir[1] - p0ToRayOrigin[1] * lineDir[0]); - // float t = CrossProduct((p0 - prd->pos), lineDir)[2] * d; - valid &= t > 1e-5f; + valid &= t > optixGetRayTmin(); float s = d * (p0ToRayOrigin[0] * prd->traceDir[1] - p0ToRayOrigin[1] * prd->traceDir[0]); - // float s = CrossProduct((p0 - prd->pos), prd->traceDir)[2] * d; valid &= s > 1e-5f && s < 1.0f - 1e-5f; if (valid) { diff --git a/gpu/pipelines/GeneralPipelineTriangle.cu b/gpu/pipelines/GeneralPipelineTriangle.cu index bf0a1ef..1bb17d2 100644 --- a/gpu/pipelines/GeneralPipelineTriangle.cu +++ b/gpu/pipelines/GeneralPipelineTriangle.cu @@ -12,8 +12,6 @@ #include "vcContext.hpp" -// #define COUNT_RAYS - using namespace viennaray::gpu; extern "C" __constant__ LaunchParams launchParams; @@ -32,8 +30,8 @@ extern "C" __global__ void __closesthit__() { // ------------- SURFACE COLLISION --------------- // unsigned callIdx; callIdx = callableIndex(launchParams.particleType, CallableSlot::COLLISION); - optixDirectCall(callIdx, sbtData, prd); + optixDirectCall(callIdx, + sbtData, prd); // ------------- REFLECTION --------------- // callIdx = callableIndex(launchParams.particleType, CallableSlot::REFLECTION); @@ -108,9 +106,5 @@ extern "C" __global__ void __raygen__() { optixReorder(hint, hintBitLength); optixInvoke(u0, u1); prd.traceDir = prd.dir; // Update traceDir for the next iteration -#ifdef COUNT_RAYS - int *counter = reinterpret_cast(launchParams.customData); - atomicAdd(counter, 1); -#endif } } diff --git a/include/viennaray/gpu/raygSource.hpp b/include/viennaray/gpu/raygSource.hpp index dd3aede..87cfcf7 100644 --- a/include/viennaray/gpu/raygSource.hpp +++ b/include/viennaray/gpu/raygSource.hpp @@ -65,16 +65,17 @@ initializeRayDirection(PerRayData &prd, const float power, __device__ __forceinline__ void initializeRayPosition(PerRayData &prd, const LaunchParams::SourcePlane &source, const uint16_t D) { - const float4 u = curand_uniform4(&prd.RNGstate); // (0,1] + const float u = curand_uniform(&prd.RNGstate); // (0,1] prd.pos[0] = - source.minPoint[0] + u.x * (source.maxPoint[0] - source.minPoint[0]); + source.minPoint[0] + u * (source.maxPoint[0] - source.minPoint[0]); if (D == 2) { prd.pos[1] = source.planeHeight; prd.pos[2] = 0.f; } else { + const float v = curand_uniform(&prd.RNGstate); // (0,1] prd.pos[1] = - source.minPoint[1] + u.y * (source.maxPoint[1] - source.minPoint[1]); + source.minPoint[1] + v * (source.maxPoint[1] - source.minPoint[1]); prd.pos[2] = source.planeHeight; } } From 67af9689b5d70d1547e06311a83fc8ef119a74eb Mon Sep 17 00:00:00 2001 From: reiter Date: Wed, 4 Feb 2026 16:36:20 +0100 Subject: [PATCH 07/10] Update triangle boundary logic --- gpu/pipelines/GeneralPipelineTriangle.cu | 23 ++++++++++++++----- include/viennaray/gpu/raygSBTRecords.hpp | 12 ++++++++-- include/viennaray/gpu/raygTraceTriangle.hpp | 6 ++--- .../viennaray/gpu/raygTriangleGeometry.hpp | 5 ++++ 4 files changed, 34 insertions(+), 12 deletions(-) diff --git a/gpu/pipelines/GeneralPipelineTriangle.cu b/gpu/pipelines/GeneralPipelineTriangle.cu index 1bb17d2..2228dca 100644 --- a/gpu/pipelines/GeneralPipelineTriangle.cu +++ b/gpu/pipelines/GeneralPipelineTriangle.cu @@ -49,12 +49,23 @@ extern "C" __global__ void __closesthit__boundary__() { prd->pos = prd->pos + prd->traceDir * optixGetRayTmax(); const unsigned int primID = optixGetPrimitiveIndex(); - unsigned dim = primID / 4; - if (launchParams.periodicBoundary) { - prd->pos[dim] = sbtData->vertex[sbtData->index[primID ^ 2][0]][dim]; - } else { - prd->dir[dim] -= 2 * prd->dir[dim]; - prd->pos[dim] = sbtData->vertex[sbtData->index[primID][0]][dim]; + // 0-3: X axis (dim 0), 4-7: Y axis (dim 1) + const unsigned int dim = primID / 4; + // 0,1,4,5 are Minimum side (0); 2,3,6,7 are Maximum side (1) + const unsigned int side = (primID & 2) >> 1; + + const int periodic = launchParams.periodicBoundary; + const float bounds[2] = {sbtData->box.minExtent[dim], + sbtData->box.maxExtent[dim]}; + + // Update Position: + // Periodic(1): opposite side (side ^ 1) + // Reflect(0): same side (side ^ 0) + prd->pos[dim] = bounds[side ^ periodic]; + + if (!launchParams.periodicBoundary) { + // Reflect direction + prd->dir[dim] = -prd->dir[dim]; } } diff --git a/include/viennaray/gpu/raygSBTRecords.hpp b/include/viennaray/gpu/raygSBTRecords.hpp index 64a0cc0..9f6b5d5 100644 --- a/include/viennaray/gpu/raygSBTRecords.hpp +++ b/include/viennaray/gpu/raygSBTRecords.hpp @@ -22,8 +22,16 @@ struct HitSBTDataDisk { struct HitSBTDataTriangle { HitSBTDataBase base; - Vec3Df *vertex; - Vec3D *index; + union { + struct { + Vec3Df *vertex; + Vec3D *index; + }; + struct { + Vec3Df minExtent; + Vec3Df maxExtent; + } box; + }; }; struct HitSBTDataLine { diff --git a/include/viennaray/gpu/raygTraceTriangle.hpp b/include/viennaray/gpu/raygTraceTriangle.hpp index 2649a90..bda8f95 100644 --- a/include/viennaray/gpu/raygTraceTriangle.hpp +++ b/include/viennaray/gpu/raygTraceTriangle.hpp @@ -91,10 +91,8 @@ template class TraceTriangle final : public Trace { HitgroupRecordTriangle boundaryHitgroupRecord = {}; optixSbtRecordPackHeader(this->boundaryHitgroupPG_, &boundaryHitgroupRecord); - boundaryHitgroupRecord.data.vertex = - (Vec3Df *)triangleGeometry_.boundaryVertexBuffer.dPointer(); - boundaryHitgroupRecord.data.index = - (Vec3D *)triangleGeometry_.boundaryIndexBuffer.dPointer(); + boundaryHitgroupRecord.data.box.minExtent = triangleGeometry_.minExtent; + boundaryHitgroupRecord.data.box.maxExtent = triangleGeometry_.maxExtent; boundaryHitgroupRecord.data.base.geometryType = 0; boundaryHitgroupRecord.data.base.isBoundary = true; diff --git a/include/viennaray/gpu/raygTriangleGeometry.hpp b/include/viennaray/gpu/raygTriangleGeometry.hpp index 5027cc4..be332c3 100644 --- a/include/viennaray/gpu/raygTriangleGeometry.hpp +++ b/include/viennaray/gpu/raygTriangleGeometry.hpp @@ -20,6 +20,9 @@ struct TriangleGeometry { CudaBuffer boundaryVertexBuffer; CudaBuffer boundaryIndexBuffer; + Vec3Df minExtent; + Vec3Df maxExtent; + // buffer that keeps the (final, compacted) accel structure CudaBuffer asBuffer; @@ -46,6 +49,8 @@ struct TriangleGeometry { mesh.maximumExtent[2] + mesh.gridDelta + sourceOffset; } launchParams.numElements = mesh.triangles.size(); + minExtent = mesh.minimumExtent; + maxExtent = mesh.maximumExtent; // 2 inputs: one for the geometry, one for the boundary std::array triangleInput{}; From 5988a908197c09ecf310bbb237f3a1fdfe79b8ec Mon Sep 17 00:00:00 2001 From: reiter Date: Wed, 11 Feb 2026 16:39:56 +0100 Subject: [PATCH 08/10] fix: direction update in CPU kernel --- include/viennaray/rayTraceKernel.hpp | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/include/viennaray/rayTraceKernel.hpp b/include/viennaray/rayTraceKernel.hpp index 7d1f8f7..31102a4 100644 --- a/include/viennaray/rayTraceKernel.hpp +++ b/include/viennaray/rayTraceKernel.hpp @@ -136,12 +136,10 @@ template class TraceKernel { pSource_->getOriginAndDirection(idx, rngState); fillRayPosition(rayHit.ray, originAndDirection[0]); if (isZero(rayDirection)) { - assert(IsNormalized(originAndDirection[1])); - fillRayDirection(rayHit.ray, originAndDirection[1]); - } else { - assert(IsNormalized(rayDirection)); - fillRayDirection(rayHit.ray, rayDirection); + rayDirection = std::move(originAndDirection[1]); } + assert(IsNormalized(rayDirection)); + fillRayDirection(rayHit.ray, rayDirection); } #ifdef VIENNARAY_USE_RAY_MASKING @@ -302,7 +300,7 @@ template class TraceKernel { } // get sticking probability (1) and reflected direction (2) - const auto stickingDirection = particle->surfaceReflection( + auto stickingDirection = particle->surfaceReflection( rayWeight, rayDirection, geomNormal, rayHit.hit.primID, geometry_.getMaterialId(rayHit.hit.primID), pGlobalData_, rngState); @@ -323,8 +321,9 @@ template class TraceKernel { } // Update ray direction and origin + rayDirection = std::move(stickingDirection.second); fillRayPosition(rayHit.ray, hitPoint); - fillRayDirection(rayHit.ray, stickingDirection.second); + fillRayDirection(rayHit.ray, rayDirection); } while (reflect); totalBoundaryHits += boundaryHits; From 405eea4e097a0da1fe5d6d3eaabaf766abc5ca51 Mon Sep 17 00:00:00 2001 From: reiter Date: Wed, 11 Feb 2026 16:57:18 +0100 Subject: [PATCH 09/10] chore: bump version --- CMakeLists.txt | 2 +- README.md | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index edf3abb..e803622 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.20 FATAL_ERROR) project( ViennaRay LANGUAGES CXX - VERSION 3.10.2) + VERSION 3.11.0) # -------------------------------------------------------------------------------------------------------- # Library switches diff --git a/README.md b/README.md index 08967b9..907815f 100644 --- a/README.md +++ b/README.md @@ -63,7 +63,7 @@ We recommend using [CPM.cmake](https://github.com/cpm-cmake/CPM.cmake) to consum * Installation with CPM ```cmake - CPMAddPackage("gh:viennatools/viennaray@3.10.2") # Use the latest release version + CPMAddPackage("gh:viennatools/viennaray@3.11.0") # Use the latest release version ``` * With a local installation From f681dd0a8aad980fbe1482abb195f1287bc03ae8 Mon Sep 17 00:00:00 2001 From: reiter Date: Wed, 11 Feb 2026 18:04:26 +0100 Subject: [PATCH 10/10] fix: tests and missing boundary hit count --- gpu/pipelines/GeneralPipelineDisk.cu | 2 + gpu/pipelines/GeneralPipelineLine.cu | 2 + gpu/pipelines/GeneralPipelineTriangle.cu | 2 + gpu/tests/boundaries/TestPipelineTriangle.cu | 66 ++++++++++++++------ gpu/tests/boundaries/boundaries.cpp | 10 +-- gpu/tests/reflections/testReflections.cu | 8 +-- tests/createRay/createRay.cpp | 14 ++--- tests/createSourceGrid/createSourceGrid.cpp | 2 +- 8 files changed, 71 insertions(+), 35 deletions(-) diff --git a/gpu/pipelines/GeneralPipelineDisk.cu b/gpu/pipelines/GeneralPipelineDisk.cu index 43b8d37..8583feb 100644 --- a/gpu/pipelines/GeneralPipelineDisk.cu +++ b/gpu/pipelines/GeneralPipelineDisk.cu @@ -145,6 +145,8 @@ extern "C" __global__ void __closesthit__boundary__() { prd->dir[axis] -= 2 * prd->dir[axis]; // reflect prd->pos[axis] = sbtData->point[primID][axis]; } + + prd->numBoundaryHits++; } extern "C" __global__ void __miss__() { getPRD()->rayWeight = 0.f; } diff --git a/gpu/pipelines/GeneralPipelineLine.cu b/gpu/pipelines/GeneralPipelineLine.cu index f21a144..51c9772 100644 --- a/gpu/pipelines/GeneralPipelineLine.cu +++ b/gpu/pipelines/GeneralPipelineLine.cu @@ -88,6 +88,8 @@ extern "C" __global__ void __closesthit__boundary__() { } else { prd->dir[0] -= 2 * prd->dir[0]; // reflect } + + prd->numBoundaryHits++; } extern "C" __global__ void __miss__() { getPRD()->rayWeight = 0.f; } diff --git a/gpu/pipelines/GeneralPipelineTriangle.cu b/gpu/pipelines/GeneralPipelineTriangle.cu index 2228dca..9f8db5d 100644 --- a/gpu/pipelines/GeneralPipelineTriangle.cu +++ b/gpu/pipelines/GeneralPipelineTriangle.cu @@ -67,6 +67,8 @@ extern "C" __global__ void __closesthit__boundary__() { // Reflect direction prd->dir[dim] = -prd->dir[dim]; } + + prd->numBoundaryHits++; } extern "C" __global__ void __miss__() { getPRD()->rayWeight = 0.f; } diff --git a/gpu/tests/boundaries/TestPipelineTriangle.cu b/gpu/tests/boundaries/TestPipelineTriangle.cu index aa6ae55..5b439f5 100644 --- a/gpu/tests/boundaries/TestPipelineTriangle.cu +++ b/gpu/tests/boundaries/TestPipelineTriangle.cu @@ -18,7 +18,7 @@ using namespace viennaray::gpu; -extern "C" __constant__ viennaray::gpu::LaunchParams launchParams; +extern "C" __constant__ LaunchParams launchParams; extern "C" __global__ void __closesthit__() { const HitSBTDataTriangle *sbtData = @@ -29,16 +29,43 @@ extern "C" __global__ void __closesthit__() { prd->tMin = optixGetRayTmax(); prd->primID = primID; - if (sbtData->base.isBoundary) { - if (launchParams.periodicBoundary) { - applyPeriodicBoundary(prd, sbtData, launchParams.D); - } else { - reflectFromBoundary(prd, sbtData, launchParams.D); - } - } else { - atomicAdd(&launchParams.resultBuffer[prd->primID], prd->rayWeight); - prd->rayWeight = 0.f; + atomicAdd(&launchParams.resultBuffer[prd->primID], prd->rayWeight); + prd->rayWeight = 0.f; + prd->load = 0.f; +} + +extern "C" __global__ void __closesthit__boundary__() { + const HitSBTDataTriangle *sbtData = + (const HitSBTDataTriangle *)optixGetSbtDataPointer(); + PerRayData *prd = getPRD(); + + // update ray position to hit point + prd->tMin = optixGetRayTmax(); + prd->pos = prd->pos + prd->dir * prd->tMin; + + const unsigned int primID = optixGetPrimitiveIndex(); + // 0-3: X axis (dim 0), 4-7: Y axis (dim 1) + const unsigned int dim = primID / 4; + // 0,1,4,5 are Minimum side (0); 2,3,6,7 are Maximum side (1) + const unsigned int side = (primID & 2) >> 1; + + const int periodic = launchParams.periodicBoundary; + const float bounds[2] = {sbtData->box.minExtent[dim], + sbtData->box.maxExtent[dim]}; + + // Update Position: + // Periodic(1): opposite side (side ^ 1) + // Reflect(0): same side (side ^ 0) + prd->pos[dim] = bounds[side ^ periodic]; + + if (!launchParams.periodicBoundary) { + // Reflect direction + prd->dir[dim] = -prd->dir[dim]; } + + prd->primID = primID; + prd->numBoundaryHits += 1; + prd->load = 1.f; } extern "C" __global__ void __miss__() { getPRD()->rayWeight = 0.f; } @@ -52,7 +79,7 @@ extern "C" __global__ void __raygen__() { // per-ray data PerRayData prd; // each ray has its own RNG state - initializeRNGState(&prd, linearLaunchIndex, launchParams.seed); + initializeRNGState(prd, linearLaunchIndex, launchParams.seed); // initialize ray position and direction if (launchParams.D == 2) { @@ -106,9 +133,9 @@ extern "C" __global__ void __raygen__() { packPointer((void *)&prd, u0, u1); while (continueRay(launchParams, prd)) { - // printf("Tracing ray %u from pos (%f, %f, %f) in dir (%f, %f, %f)\n", - // linearLaunchIndex, prd.pos[0], prd.pos[1], prd.pos[2], prd.dir[0], - // prd.dir[1], prd.dir[2]); + printf("Tracing ray %u from pos (%f, %f, %f) in dir (%f, %f, %f)\n", + linearLaunchIndex, prd.pos[0], prd.pos[1], prd.pos[2], prd.dir[0], + prd.dir[1], prd.dir[2]); optixTrace(launchParams.traversable, // traversable GAS make_float3(prd.pos[0], prd.pos[1], prd.pos[2]), // origin make_float3(prd.dir[0], prd.dir[1], prd.dir[2]), // direction @@ -121,9 +148,12 @@ extern "C" __global__ void __raygen__() { 1, // SBT stride 0, // missSBTIndex u0, u1); // Payload - // printf("Ray %u hit primID %u with tMin %f\n", linearLaunchIndex, - // prd.primID, - // prd.tMin); - prd.numReflections++; + if (prd.load > 0.f) { + printf("Ray %u hit BOUNDARY primID %u with tMin %f\n", linearLaunchIndex, + prd.primID, prd.tMin); + } else { + printf("Ray %u hit SURFACE primID %u with tMin %f\n", linearLaunchIndex, + prd.primID, prd.tMin); + } } } diff --git a/gpu/tests/boundaries/boundaries.cpp b/gpu/tests/boundaries/boundaries.cpp index 3194186..737b453 100644 --- a/gpu/tests/boundaries/boundaries.cpp +++ b/gpu/tests/boundaries/boundaries.cpp @@ -68,11 +68,11 @@ int main() { tracer.apply(); auto flux = tracer.getFlux(0, 0, 0); - // std::cout << "Flux values at each triangle:" << std::endl; - // for (size_t i = 0; i < flux.size(); ++i) { - // std::cout << "Triangle " << i << ": " << flux[i] << std::endl; - // } - // rayInternal::writeVTP(mesh, "boundary_3d_reflective.vtp", flux); + std::cout << "Flux values at each triangle:" << std::endl; + for (size_t i = 0; i < flux.size(); ++i) { + std::cout << "Triangle " << i << ": " << flux[i] << std::endl; + } + rayInternal::writeVTP(mesh, "boundary_3d_reflective.vtp", flux); VC_TEST_ASSERT(flux[3] > 0.f); VC_TEST_ASSERT(flux[5] > 0.f); diff --git a/gpu/tests/reflections/testReflections.cu b/gpu/tests/reflections/testReflections.cu index 7198e27..a7bc0a5 100644 --- a/gpu/tests/reflections/testReflections.cu +++ b/gpu/tests/reflections/testReflections.cu @@ -16,10 +16,10 @@ extern "C" __global__ void test_diffuse(viennacore::Vec3Df inDir, for (; tidx < numResults; tidx += stride) { viennaray::gpu::PerRayData prd; - initializeRNGState(&prd, tidx, 0); + initializeRNGState(prd, tidx, 0); prd.dir = inDir; - diffuseReflection(&prd, normal, 3); + diffuseReflection(&prd, normal); results[tidx] = prd.dir; } } @@ -35,10 +35,10 @@ extern "C" __global__ void test_coned_cosine(viennacore::Vec3Df inDir, for (; tidx < numResults; tidx += stride) { viennaray::gpu::PerRayData prd; - initializeRNGState(&prd, tidx, 0); + initializeRNGState(prd, tidx, 0); prd.dir = inDir; - conedCosineReflection(&prd, normal, coneAngle, 3); + conedCosineReflection(&prd, normal, coneAngle); results[tidx] = prd.dir; } } diff --git a/tests/createRay/createRay.cpp b/tests/createRay/createRay.cpp index a62dc82..17090a8 100644 --- a/tests/createRay/createRay.cpp +++ b/tests/createRay/createRay.cpp @@ -51,7 +51,7 @@ int main() { for (size_t i = 0; i < 10; ++i) { auto originAndDirection = source.getOriginAndDirection(0, rngState); rayInternal::fillRayPosition(rayHit.ray, originAndDirection[0]); - rayInternal::fillRayDirection(rayHit.ray, originAndDirection[1]); + rayInternal::fillRayDirection(rayHit.ray, originAndDirection[1]); VC_TEST_ASSERT(rayHit.ray.dir_z < 0.) VC_TEST_ASSERT_ISCLOSE(rayHit.ray.org_z, (1. + 2 * gridDelta), eps) } @@ -73,7 +73,7 @@ int main() { for (size_t i = 0; i < 10; ++i) { auto originAndDirection = source.getOriginAndDirection(0, rngState); rayInternal::fillRayPosition(rayHit.ray, originAndDirection[0]); - rayInternal::fillRayDirection(rayHit.ray, originAndDirection[1]); + rayInternal::fillRayDirection(rayHit.ray, originAndDirection[1]); VC_TEST_ASSERT(rayHit.ray.dir_z > 0.) VC_TEST_ASSERT_ISCLOSE(rayHit.ray.org_z, (-1. - 2 * gridDelta), eps) } @@ -95,7 +95,7 @@ int main() { for (size_t i = 0; i < 10; ++i) { auto originAndDirection = source.getOriginAndDirection(0, rngState); rayInternal::fillRayPosition(rayHit.ray, originAndDirection[0]); - rayInternal::fillRayDirection(rayHit.ray, originAndDirection[1]); + rayInternal::fillRayDirection(rayHit.ray, originAndDirection[1]); VC_TEST_ASSERT(rayHit.ray.dir_x < 0.) VC_TEST_ASSERT_ISCLOSE(rayHit.ray.org_x, (1. + 2 * gridDelta), eps) } @@ -117,7 +117,7 @@ int main() { for (size_t i = 0; i < 10; ++i) { auto originAndDirection = source.getOriginAndDirection(0, rngState); rayInternal::fillRayPosition(rayHit.ray, originAndDirection[0]); - rayInternal::fillRayDirection(rayHit.ray, originAndDirection[1]); + rayInternal::fillRayDirection(rayHit.ray, originAndDirection[1]); VC_TEST_ASSERT(rayHit.ray.dir_x > 0.) VC_TEST_ASSERT_ISCLOSE(rayHit.ray.org_x, (-1. - 2 * gridDelta), eps) } @@ -139,7 +139,7 @@ int main() { for (size_t i = 0; i < 10; ++i) { auto originAndDirection = source.getOriginAndDirection(0, rngState); rayInternal::fillRayPosition(rayHit.ray, originAndDirection[0]); - rayInternal::fillRayDirection(rayHit.ray, originAndDirection[1]); + rayInternal::fillRayDirection(rayHit.ray, originAndDirection[1]); VC_TEST_ASSERT(rayHit.ray.dir_y < 0.) VC_TEST_ASSERT_ISCLOSE(rayHit.ray.org_y, (1. + 2 * gridDelta), eps) } @@ -161,7 +161,7 @@ int main() { for (size_t i = 0; i < 10; ++i) { auto originAndDirection = source.getOriginAndDirection(0, rngState); rayInternal::fillRayPosition(rayHit.ray, originAndDirection[0]); - rayInternal::fillRayDirection(rayHit.ray, originAndDirection[1]); + rayInternal::fillRayDirection(rayHit.ray, originAndDirection[1]); VC_TEST_ASSERT(rayHit.ray.dir_y > 0.) VC_TEST_ASSERT_ISCLOSE(rayHit.ray.org_y, (-1. - 2 * gridDelta), eps) } @@ -186,7 +186,7 @@ int main() { for (size_t i = 0; i < 10; ++i) { auto originAndDirection = source.getOriginAndDirection(0, rngState); rayInternal::fillRayPosition(rayHit.ray, originAndDirection[0]); - rayInternal::fillRayDirection(rayHit.ray, originAndDirection[1]); + rayInternal::fillRayDirection(rayHit.ray, originAndDirection[1]); VC_TEST_ASSERT(rayHit.ray.dir_z < 0.) VC_TEST_ASSERT_ISCLOSE(rayHit.ray.org_z, (1. + 2 * gridDelta), eps) } diff --git a/tests/createSourceGrid/createSourceGrid.cpp b/tests/createSourceGrid/createSourceGrid.cpp index a785930..9b51270 100644 --- a/tests/createSourceGrid/createSourceGrid.cpp +++ b/tests/createSourceGrid/createSourceGrid.cpp @@ -38,7 +38,7 @@ int main() { for (size_t i = 0; i < numGridPoints; ++i) { auto originAndDirection = source.getOriginAndDirection(i, rngState); rayInternal::fillRayPosition(rayHit.ray, originAndDirection[0]); - rayInternal::fillRayDirection(rayHit.ray, originAndDirection[1]); + rayInternal::fillRayDirection(rayHit.ray, originAndDirection[1]); VC_TEST_ASSERT(rayHit.ray.dir_z < 0.) VC_TEST_ASSERT_ISCLOSE(rayHit.ray.org_z, (1. + 2 * gridDelta), eps)