Update pipelines (#79)

* Enable custom direction basis and optimize data types

* Add ray load

* Remove TIndex array from PRD

* Correct numBoundaryHits

* Add option to precompute all normals

* Fix type

* Small fixes

* Bump version

* Remove unnecessary enums

Author: tobre1

CMakeLists.txt

@@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.20 FATAL_ERROR)
 project(
   ViennaRay
   LANGUAGES CXX
-  VERSION 3.6.0)
+  VERSION 3.6.1)
 
 # --------------------------------------------------------------------------------------------------------
 # Library switches

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.6.0")
+  CPMAddPackage("gh:viennatools/viennaray@3.6.1")
   ```
 
 * With a local installation

gpu/examples/CMakeLists.txt

@@ -7,7 +7,8 @@ configure_file(Resources/trenchMesh.dat ${CMAKE_CURRENT_BINARY_DIR}/trenchMesh.d
 add_dependencies(ViennaRay-GPU_Examples ${target_name})
 
 add_gpu_executable(GPU_trenchDisks target_name trenchDisks.cpp)
-configure_file(${CMAKE_SOURCE_DIR}/examples/trench/trenchGrid3D.dat ${CMAKE_CURRENT_BINARY_DIR}/trenchGrid3D.dat COPYONLY)
+configure_file(${CMAKE_SOURCE_DIR}/examples/trench/trenchGrid3D.dat
+               ${CMAKE_CURRENT_BINARY_DIR}/trenchGrid3D.dat COPYONLY)
 add_dependencies(ViennaRay-GPU_Examples ${target_name})
 
 add_gpu_executable(GPU_trenchLines target_name trenchLines.cpp)

gpu/include/raygBoundary.hpp

@@ -16,6 +16,7 @@ __device__ __inline__ void reflectFromBoundary(viennaray::gpu::PerRayData *prd,
                                                const int D) {
   using namespace viennacore;
   const unsigned int primID = optixGetPrimitiveIndex();
+  prd->numBoundaryHits++;
 
   if constexpr (std::is_same<SBTData, viennaray::gpu::HitSBTDataDisk>::value) {
     prd->pos =
@@ -32,7 +33,6 @@ __device__ __inline__ void reflectFromBoundary(viennaray::gpu::PerRayData *prd,
     unsigned dim = primID / 4;
     prd->dir[dim] -= 2 * prd->dir[dim];
     prd->pos[dim] = hsd->vertex[hsd->index[primID][0]][dim];
-    prd->numBoundaryHits++;
   } else if constexpr (std::is_same<SBTData,
                                     viennaray::gpu::HitSBTDataLine>::value) {
     prd->pos = prd->pos + prd->dir * (optixGetRayTmax());
@@ -47,6 +47,7 @@ applyPeriodicBoundary(viennaray::gpu::PerRayData *prd, const SBTData *hsd,
                       const int D) {
   using namespace viennacore;
   const unsigned int primID = optixGetPrimitiveIndex();
+  prd->numBoundaryHits++;
 
   if constexpr (std::is_same<SBTData, viennaray::gpu::HitSBTDataDisk>::value) {
     prd->pos =
@@ -66,7 +67,6 @@ applyPeriodicBoundary(viennaray::gpu::PerRayData *prd, const SBTData *hsd,
     prd->pos = prd->pos + prd->dir * optixGetRayTmax();
     unsigned dim = primID / 4;
     prd->pos[dim] = hsd->vertex[hsd->index[primID ^ 2][0]][dim];
-    prd->numBoundaryHits++;
   } else if constexpr (std::is_same<SBTData,
                                     viennaray::gpu::HitSBTDataLine>::value) {
     prd->pos = prd->pos + prd->dir * (optixGetRayTmax());

gpu/include/raygLaunchParams.hpp

@@ -15,32 +15,34 @@ __both__ __forceinline__ unsigned callableIndex(unsigned p, CallableSlot s) {
 
 struct LaunchParams {
   float *resultBuffer;
+
   float rayWeightThreshold = 0.1f;
+  float tThreshold = 0.5f;
+
   unsigned int seed = 0;
-  unsigned int numElements;
-  unsigned int *dataPerParticle; // to determine result buffer index
   bool periodicBoundary = false;
-  unsigned int maxBoundaryHits = 100;
-  unsigned int particleIdx = 0;
-  unsigned particleType = 0;
-  float gridDelta = 1.f;
-  float tThreshold = 0.5f;
 
-  int D = 3; // Dimension
+  unsigned int numElements;      // to determine result buffer index
+  unsigned int *dataPerParticle; // to determine result buffer index
+
+  unsigned maxBoundaryHits = 100;
+  uint8_t particleIdx = 0;
+  uint8_t particleType = 0;
+  uint8_t D = 3; // Dimension
 
-  // std::unordered_map<int, float> sticking;
-  int *materialIds;
-  float *materialSticking;
   float sticking = 1.f;
   float cosineExponent = 1.f;
+  int *materialIds;
+  float *materialSticking;
   void *customData;
 
   // source plane params
-  struct {
+  struct SourcePlane {
     viennacore::Vec2Df minPoint;
     viennacore::Vec2Df maxPoint;
     float planeHeight;
     std::array<viennacore::Vec3Df, 3> directionBasis;
+    bool customDirectionBasis = false;
   } source;
 
   OptixTraversableHandle traversable;

gpu/include/raygLineGeometry.hpp

@@ -15,6 +15,7 @@ template <typename NumericType, int D = 3> struct LineGeometry {
   // geometry
   CudaBuffer geometryNodesBuffer;
   CudaBuffer geometryLinesBuffer;
+  CudaBuffer geometryNormalsBuffer;
 
   // boundary
   CudaBuffer boundaryNodesBuffer;
@@ -52,6 +53,7 @@ template <typename NumericType, int D = 3> struct LineGeometry {
     // upload the model to the device: the builder
     geometryNodesBuffer.allocUpload(mesh.nodes);
     geometryLinesBuffer.allocUpload(mesh.lines);
+    geometryNormalsBuffer.allocUpload(mesh.normals);
 
     // create local variables, because we need a *pointer* to the
     // device pointers

gpu/include/raygMesh.hpp

@@ -13,6 +13,7 @@ using namespace viennacore;
 struct LineMesh {
   std::vector<Vec3Df> nodes;
   std::vector<Vec2D<unsigned>> lines;
+  std::vector<Vec3Df> normals;
 
   Vec3Df minimumExtent;
   Vec3Df maximumExtent;
@@ -22,6 +23,7 @@ struct LineMesh {
 struct TriangleMesh {
   std::vector<Vec3Df> nodes;
   std::vector<Vec3D<unsigned>> triangles;
+  std::vector<Vec3Df> normals;
 
   Vec3Df minimumExtent;
   Vec3Df maximumExtent;

gpu/include/raygPerRayData.hpp

@@ -12,25 +12,31 @@
 
 namespace viennaray::gpu {
 
+// Per-ray data structure associated with each ray. Should be kept small to
+// optimize memory usage and performance.
 struct PerRayData {
-  float rayWeight = 1.f;
+  // Position and direction
   viennacore::Vec3Df pos;
   viennacore::Vec3Df dir;
 
+  // Simulation specific data
+  float rayWeight = 1.f;
+  float energy = 0.f;
+  float load = 0.f;
+
+  // RNG
   RNGState RNGstate;
 
-  float energy = 0.f;
+  // Hit data
   unsigned int numBoundaryHits = 0;
-
-  unsigned primID = 0;
-  float tMin = 1e20f;
+  unsigned int primID = 0; // primID of closest hit
+  float tMin = 1e20f;      // distance to closest hit
 
   // Variables for neighbor intersections (overlapping disks and lines)
-  int TIndex[MAX_NEIGHBORS];    // Indices of neighbor hits
-  int ISCount = 0;              // Number of hits starting from 1
-  int tempCount = 0;            // total intersections recorded
-  float tValues[MAX_NEIGHBORS]; // all intersection distances
-  int primIDs[MAX_NEIGHBORS];   // their primitive IDs
+  uint8_t ISCount = 0;                 // Number of hits starting from 1
+  uint8_t totalCount = 0;              // total intersections recorded
+  float tValues[MAX_NEIGHBORS];        // all intersection distances
+  unsigned int primIDs[MAX_NEIGHBORS]; // their primitive IDs
   bool hitFromBack = false;
 };
 

gpu/include/raygRNG.hpp

@@ -37,10 +37,9 @@ __device__ __inline__ float getNextRand(viennaray::gpu::RNGState *state) {
 }
 
 __device__ __inline__ float getNormalDistRand(viennaray::gpu::RNGState *state) {
-  float u0 = curand_uniform(state);
-  float u1 = curand_uniform(state);
-  float r = sqrtf(-2.f * logf(u0));
-  float theta = 2.f * M_PIf * u1;
+  float4 u0 = curand_uniform4(state);
+  float r = sqrtf(-2.f * logf(u0.x));
+  float theta = 2.f * M_PIf * u0.y;
   return r * sinf(theta);
 }
 #endif

gpu/include/raygReflection.hpp

@@ -13,11 +13,12 @@ using namespace viennaray::gpu;
 #ifdef __CUDACC__
 __device__ __inline__ viennacore::Vec3Df
 computeNormal(const void *sbtData, const unsigned int primID) {
+  using namespace viennacore;
   const HitSBTDataBase *baseData =
       reinterpret_cast<const HitSBTDataBase *>(sbtData);
   switch (baseData->geometryType) {
   case 0: {
-    using namespace viennacore;
+    // Triangles
     const HitSBTDataTriangle *sbt =
         reinterpret_cast<const HitSBTDataTriangle *>(sbtData);
     const Vec3D<unsigned> &index = sbt->index[primID];
@@ -27,12 +28,11 @@ computeNormal(const void *sbtData, const unsigned int primID) {
     return Normalize<float, 3>(CrossProduct<float>(B - A, C - A));
   } break;
   case 1: {
-    const HitSBTDataDisk *sbt =
-        reinterpret_cast<const HitSBTDataDisk *>(sbtData);
-    return sbt->normal[primID];
+    // Disks
+    return baseData->normal[primID];
   } break;
   case 2: {
-    using namespace viennacore;
+    // Lines
     const HitSBTDataLine *sbt =
         reinterpret_cast<const HitSBTDataLine *>(sbtData);
     Vec3Df p0 = sbt->nodes[sbt->lines[primID][0]];
@@ -48,6 +48,11 @@ computeNormal(const void *sbtData, const unsigned int primID) {
   }
 }
 
+__device__ __forceinline__ viennacore::Vec3Df
+getNormal(const void *sbtData, const unsigned int primID) {
+  return reinterpret_cast<const HitSBTDataBase *>(sbtData)->normal[primID];
+}
+
 static __device__ __forceinline__ void
 specularReflection(viennaray::gpu::PerRayData *prd,
                    const viennacore::Vec3Df &geoNormal) {
@@ -163,7 +168,7 @@ PickRandomPointOnUnitSphere(viennaray::gpu::RNGState *state) {
 
 static __device__ void diffuseReflection(viennaray::gpu::PerRayData *prd,
                                          const viennacore::Vec3Df &geoNormal,
-                                         const int D) {
+                                         const uint8_t D) {
   using namespace viennacore;
 #ifndef VIENNARAY_TEST
   prd->pos = prd->pos + prd->tMin * prd->dir;
@@ -178,7 +183,7 @@ static __device__ void diffuseReflection(viennaray::gpu::PerRayData *prd,
 }
 
 static __device__ void diffuseReflection(viennaray::gpu::PerRayData *prd,
-                                         const int D) {
+                                         const uint8_t D) {
   using namespace viennacore;
 
   const viennaray::gpu::HitSBTDataDisk *sbtData =
@@ -190,7 +195,7 @@ static __device__ void diffuseReflection(viennaray::gpu::PerRayData *prd,
 static __device__ __forceinline__ void
 conedCosineReflection(viennaray::gpu::PerRayData *prd,
                       const viennacore::Vec3Df &geomNormal,
-                      const float maxConeAngle, const int D) {
+                      const float maxConeAngle, const uint8_t D) {
   using namespace viennacore;
 
   // TODO: Is this needed? Done in the CPU version