Fix periodic boundary conditions

CMakeLists.txt

@@ -15,7 +15,7 @@ option(VIENNAPS_BUILD_PYTHON "Build python bindings" OFF)
 option(VIENNAPS_PACKAGE_PYTHON "Build python bindings with intent to publish wheel" OFF)
 
 option(VIENNAPS_BUILD_APPLICATION "Build application" OFF)
-option(VIENNAPS_BUILD_GPU_APPLICATION "Build ViennaPS GPU application." OFF)
+option(VIENNAPS_BUILD_GPU_APPLICATION "Build ViennaPS GPU application." ON)
 option(VIENNAPS_USE_DOUBLE "Use double precision floating point numbers in application" OFF)
 
 # --------------------------------------------------------------------------------------------------------

gpu/CMakeLists.txt

@@ -15,7 +15,7 @@ set_property(DIRECTORY PROPERTY EXCLUDE_FROM_ALL TRUE)
 #       $ENV{OptiX_INSTALL_DIR}
 #       CACHE PATH "Path to OptiX installation.")
 # endif()
-set(OptiX_INSTALL_DIR "/home/reiter/OptiX" CACHE PATH "Path to OptiX installed location.")
+set(OptiX_INSTALL_DIR "/home/tobias/OptiX" CACHE PATH "Path to OptiX installed location.")
 set(OptiX_INCLUDE "${OptiX_INSTALL_DIR}/include")
 
 include(${CMAKE_SOURCE_DIR}/cmake/configure_optix.cmake)

gpu/pipelines/FaradayCagePipeline.cu

@@ -27,7 +27,6 @@ extern "C" __global__ void __closesthit__ion()
 {
   const HitSBTData *sbtData = (const HitSBTData *)optixGetSbtDataPointer();
   PerRayData *prd = (PerRayData *)getPRD<PerRayData>();
-  // printf("Direction: %f %f %f\n", prd->dir[0], prd->dir[1], prd->dir[2]);
 
   if (sbtData->isBoundary)
   {
@@ -58,7 +57,7 @@ extern "C" __global__ void __closesthit__ion()
     // ---------- REFLECTION ------------ //
     // prd->rayWeight -= prd->rayWeight * params.sticking;
     // conedCosineReflection(prd, (float)(M_PIf / 2.f - min(incAngle, minAngle)), geomNormal);
-    // specularReflection(prd);
+    specularReflection(prd);
   }
 }
 

gpu/rayTracing/curtBoundary.hpp

@@ -9,56 +9,56 @@
 // this can only get compiled if included in a cuda kernel
 #ifdef __CUDACC__
 __device__ __inline__ viennacore::Vec3Df
-computeNormal(const viennaps::gpu::HitSBTData *sbt, const unsigned int primID)
-{
+computeNormal(const viennaps::gpu::HitSBTData *sbt, const unsigned int primID) {
   using namespace viennacore;
-  const Vec3D<int> index = sbt->index[primID];
+  const Vec3D<unsigned> index = sbt->index[primID];
   const Vec3Df &A = sbt->vertex[index[0]];
   const Vec3Df &B = sbt->vertex[index[1]];
   const Vec3Df &C = sbt->vertex[index[2]];
   return Normalize(CrossProduct(B - A, C - A));
 }
 
-__device__ __inline__ void reflectFromBoundary(viennaps::gpu::PerRayData *prd)
-{
+__device__ __inline__ void reflectFromBoundary(viennaps::gpu::PerRayData *prd) {
   using namespace viennacore;
   const unsigned int primID = optixGetPrimitiveIndex();
   prd->pos = prd->pos + prd->dir * optixGetRayTmax();
 
-  if (primID < 4)
-  {
+  if (primID < 4) {
     prd->dir[0] -= 2 * prd->dir[0];
-  }
-  else
-  {
+  } else {
     prd->dir[1] -= 2 * prd->dir[1];
   }
 }
 
 __device__ __inline__ void
 applyPeriodicBoundary(viennaps::gpu::PerRayData *prd,
-                      const viennaps::gpu::HitSBTData *hsd)
-{
+                      const viennaps::gpu::HitSBTData *hsd) {
+  using namespace viennacore;
   const unsigned int primID = optixGetPrimitiveIndex();
+  prd->pos = prd->pos + optixGetRayTmax() * prd->dir;
 
-  // printf("Hit %u at [%f, %f, %f].", primID, prd->pos[0], prd->pos[1],
-  //        prd->pos[2]);
-  if (primID == 0 || primID == 1)
-  {
-    prd->pos[0] = hsd->vertex[hsd->index[0][0]][0]; // set to x min
-  }
-  else if (primID == 2 || primID == 3)
+  // const uint3 idx = optixGetLaunchIndex();
+  // const uint3 dims = optixGetLaunchDimensions();
+  // const unsigned int linearLaunchIndex =
+  //     idx.x + idx.y * dims.x + idx.z * dims.x * dims.y;
+
+  // printf("Ray %u Hit %u at [%f, %f, %f].\n", linearLaunchIndex, primID,
+  //        prd->pos[0], prd->pos[1], prd->pos[2]);
+
+  if (primID == 0 || primID == 1) // x min
   {
     prd->pos[0] = hsd->vertex[hsd->index[2][0]][0]; // set to x max
-  }
-  else if (primID == 4 || primID == 5)
+  } else if (primID == 2 || primID == 3)            // x max
   {
-    prd->pos[1] = hsd->vertex[hsd->index[4][0]][1]; // set to y min
-  }
-  else if (primID == 6 || primID == 7)
+    prd->pos[0] = hsd->vertex[hsd->index[0][0]][0]; // set to x min
+  } else if (primID == 4 || primID == 5)            // y min
   {
     prd->pos[1] = hsd->vertex[hsd->index[6][0]][1]; // set to y max
+  } else if (primID == 6 || primID == 7)            // y max
+  {
+    prd->pos[1] = hsd->vertex[hsd->index[4][0]][1]; // set to y min
   }
-  // printf("  Now at [%f, %f, %f].\n", prd->pos[0], prd->pos[1], prd->pos[2]);
+  // printf("Ray %u Now at [%f, %f, %f].\n", linearLaunchIndex, prd->pos[0],
+  //        prd->pos[1], prd->pos[2]);
 }
 #endif

gpu/rayTracing/curtGeometry.hpp

@@ -76,16 +76,14 @@ template <typename T, int D> struct Geometry {
     CUdeviceptr d_geoIndices = geometryIndexBuffer.d_pointer();
 
     triangleInput[0].triangleArray.vertexFormat = OPTIX_VERTEX_FORMAT_FLOAT3;
-    triangleInput[0].triangleArray.vertexStrideInBytes =
-        sizeof(std::array<float, D>);
+    triangleInput[0].triangleArray.vertexStrideInBytes = sizeof(Vec3Df);
     triangleInput[0].triangleArray.numVertices =
         (unsigned int)mesh->nodes.size();
     triangleInput[0].triangleArray.vertexBuffers = &d_geoVertices;
 
     triangleInput[0].triangleArray.indexFormat =
         OPTIX_INDICES_FORMAT_UNSIGNED_INT3;
-    triangleInput[0].triangleArray.indexStrideInBytes =
-        sizeof(std::array<unsigned, 3>);
+    triangleInput[0].triangleArray.indexStrideInBytes = sizeof(Vec3D<unsigned>);
     triangleInput[0].triangleArray.numIndexTriplets =
         (unsigned int)mesh->triangles.size();
     triangleInput[0].triangleArray.indexBuffer = d_geoIndices;
@@ -121,7 +119,7 @@ template <typename T, int D> struct Geometry {
 
     triangleInput[1].triangleArray.indexFormat =
         OPTIX_INDICES_FORMAT_UNSIGNED_INT3;
-    triangleInput[1].triangleArray.indexStrideInBytes = sizeof(Vec3D<int>);
+    triangleInput[1].triangleArray.indexStrideInBytes = sizeof(Vec3D<unsigned>);
     triangleInput[1].triangleArray.numIndexTriplets =
         (int)boundaryMesh.index.size();
     triangleInput[1].triangleArray.indexBuffer = d_boundIndices;
@@ -216,15 +214,15 @@ template <typename T, int D> struct Geometry {
     boundaryMesh.vertex.push_back(Vec3Df{bbMin[0], bbMax[1], bbMax[2]});
 
     // x min max
-    boundaryMesh.index.push_back(Vec3D<int>{0, 3, 7});
-    boundaryMesh.index.push_back(Vec3D<int>{0, 7, 4});
-    boundaryMesh.index.push_back(Vec3D<int>{6, 2, 1});
-    boundaryMesh.index.push_back(Vec3D<int>{6, 1, 5});
+    boundaryMesh.index.push_back(Vec3D<int>{0, 3, 7}); // 0
+    boundaryMesh.index.push_back(Vec3D<int>{0, 7, 4}); // 1
+    boundaryMesh.index.push_back(Vec3D<int>{6, 2, 1}); // 2
+    boundaryMesh.index.push_back(Vec3D<int>{6, 1, 5}); // 3
     // y min max
-    boundaryMesh.index.push_back(Vec3D<int>{0, 4, 5});
-    boundaryMesh.index.push_back(Vec3D<int>{0, 5, 1});
-    boundaryMesh.index.push_back(Vec3D<int>{6, 7, 3});
-    boundaryMesh.index.push_back(Vec3D<int>{6, 3, 2});
+    boundaryMesh.index.push_back(Vec3D<int>{0, 4, 5}); // 4
+    boundaryMesh.index.push_back(Vec3D<int>{0, 5, 1}); // 5
+    boundaryMesh.index.push_back(Vec3D<int>{6, 7, 3}); // 6
+    boundaryMesh.index.push_back(Vec3D<int>{6, 3, 2}); // 7
 
     boundaryMesh.minCoords = bbMin;
     boundaryMesh.maxCoords = bbMax;

gpu/rayTracing/curtSBTRecords.hpp

@@ -9,14 +9,14 @@ namespace gpu {
 
 struct HitSBTData {
   viennacore::Vec3Df *vertex;
-  viennacore::Vec3D<int> *index;
+  viennacore::Vec3D<unsigned> *index;
   bool isBoundary;
   void *cellData;
 };
 
 struct HitSBTDiskData {
   viennacore::Vec3Df *vertex;
-  viennacore::Vec3D<int> *index;
+  viennacore::Vec3D<unsigned> *index;
   bool isBoundary;
   void *cellData;
 };

gpu/rayTracing/curtTracer.hpp

@@ -582,7 +582,7 @@ template <class T, int D> class Tracer {
     geometryHitgroupRecord.data.vertex =
         (Vec3Df *)geometry.geometryVertexBuffer.d_pointer();
     geometryHitgroupRecord.data.index =
-        (Vec3D<int> *)geometry.geometryIndexBuffer.d_pointer();
+        (Vec3D<unsigned> *)geometry.geometryIndexBuffer.d_pointer();
     geometryHitgroupRecord.data.isBoundary = false;
     geometryHitgroupRecord.data.cellData = (void *)cellDataBuffer.d_pointer();
     hitgroupRecords.push_back(geometryHitgroupRecord);
@@ -593,7 +593,7 @@ template <class T, int D> class Tracer {
     boundaryHitgroupRecord.data.vertex =
         (Vec3Df *)geometry.boundaryVertexBuffer.d_pointer();
     boundaryHitgroupRecord.data.index =
-        (Vec3D<int> *)geometry.boundaryIndexBuffer.d_pointer();
+        (Vec3D<unsigned> *)geometry.boundaryIndexBuffer.d_pointer();
     boundaryHitgroupRecord.data.isBoundary = true;
     hitgroupRecords.push_back(boundaryHitgroupRecord);